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Summary 

Models that predict the failure probability of brittle materi- 
als under multiaxial loading have been developed by authors 
such as Batdorf, Evans, and Matsuo. These “unit-sphere” 
models assume that the strength-controlling flaws are ran- 
domly oriented, noninteracting planar microcracks of specified 
geometry but of variable size. This methodology has been 
extended to predict the multiaxial stochastic strength response 
of transversely isotropic brittle materials, including polymer 
matrix composites (PMCs), by considering (1) flaw- 
orientation anisotropy, whereby a preexisting microcrack has a 
higher likelihood of being oriented in one direction over 
another direction, and (2) critical strength, or K ic orientation 
anisotropy, whereby the level of critical strength or fracture 
toughness for mode I crack propagation K lr changes with 
regard to the orientation of the microstructure. In this report, 
results from finite element analysis of a fiber-reinforced- 
matrix unit cell were used with the unit-sphere model to 
predict the biaxial strength response of a unidirectional PMC 
previously reported from the World-Wide Failure Exercise. 
Results for nuclear-grade graphite materials under biaxial 
loading are also shown for comparison. This effort was 
successful in predicting the multiaxial strength response for 
the chosen problems. Findings regarding stress-state interac- 
tions and failure modes also are provided. 

1.0 Introduction 

The main purpose of this report is to investigate how “unit- 
sphere” methodology can be used as a failure criterion inside 
finite-element- or micromechanics-based software codes for 
individual brittle (or quasi-brittle) composite-material constit- 
uents to predict the overall strength response. The term “unit- 
sphere” as used herein refers to the models that were devel- 
oped by Batdorf and Crose (1974), Batdorf and Fleinisch 
(1978), Evans (1978), and Matsuo (1981) to predict the 
probability of failure of brittle materials under multiaxial 
loading. These models use a unit radius sphere representing 
the random orientation of flaws to calculate the effect of 
multiaxial stresses on material reliability. This approach 
assumes that the strength-controlling flaws are randomly 
oriented, noninteracting planar microcracks of specified 
geometry but of variable size. Fracture mechanics relation- 


ships for mixed-mode (modes I, II, and III) crack growth, 
combined with the weakest-link theory and integration over 
the surface area of a unit radius sphere representing all 
possible orientations of microcracks, are used to calculate the 
material probability of failure. This unit-sphere methodology 
was originally introduced within a statistical theory of brittle 
material strength by Weibull (1939), though without consider- 
ation for the mechanics of crack growth. 

This report describes the development of the unit-sphere 
methodology for generalized anisotropic strength response 
(in this case, for transverse isotropy) and the use of the 
methodology to predict the biaxial strength response of a 
unidirectional polymer matrix composite (PMC) previously 
reported from the World-Wide Failure Exercise (WWFE) from 
Flinton et al. (2004). Also shown are results for nuclear-grade 
graphite materials under biaxial loading. These two examples 
represent an application to a homogeneous material (i.e., 
graphite) and a heterogeneous material system (i.e., a fiber- 
reinforced composite) under multiaxial loading. A follow-on 
phase of this effort, not described here, allows this technology 
to work with NASA’s MAC/GMC micromechanics analysis 
code (MAC) of Bednarcyk and Arnold (2002), which is based 
on the generalized method of cells (GMC) family of 
micromechanics theories, including doubly and triply periodic 
versions of the GMC (Aboudi (1995)) and the High-F idelity 
Generalized Method of Cells (F1FGMC) (Aboudi et al. (2003)). 
This incorporation will allow the full exercise of the unit-sphere 
methodology, including incremental time/load steps and fatigue 
analysis (as described in Nemeth et al. (2005)), to predict the 
durability (strength and lifetime) of composite laminates and 
woven composite structures. 

Central to this report is the assumption that a composite 
material is controlled by various independent failure modes 
and that these failure modes follow a Weibull distribution. The 
product of the probability of survival of the various failure 
modes yields the probability of survival (the reliability) of the 
composite. This is a weakest-link assumption. The weakest- 
link theory (WLT) and the Weibull distribution predict a size 
effect (whereby larger components have a lower average 
strength than equally loaded smaller components). Size-effect 
studies are often used by researchers as a proxy to test WLT 
behavior and the applicability of the Weibull strength theory 
for a material. 

The review article by Wisnom (1999) examines the then- 
available experimental data for PMCs regarding size effect 
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and Weibull distribution for various composite failure modes. 
Various defect populations and factors that influence the size 
effect are discussed, as are the strengths and weaknesses of the 
Weibull approach. Wisnom (1999) indicated that there was 
broad evidence from the literature for size effects for all 
failure modes, although there were also inconsistencies and 
complications in the experimental record. Generally Wisnom 
(1999) concluded that the size effect was consistent with a 
Weibull modulus between 13 and 29 and that size effect was 
larger for matrix-dominated failure, for which the experi- 
mental data correlated reasonably well with the trend pre- 
dicted by the Weibull strength theory. For ceramic matrix 
composites (CMCs), size effect has been studied by others — 
for example, Calard and Lamon (2002), Lamon (2001), and 
Nozawa et al. (2002), although the experimental record is 
relatively meager in comparison to that for PMCs. These 
papers are just a few noteworthy examples from the literature 
on the subject of size effect and strength scatter in composites. 

It is clear that the phenomenon of damage accumulation and 
interaction, as opposed to standard WLT, complicates the 
situation regarding size effect. However, by focusing instead on 
the failure behavior of the individual composite material 
constituents of fiber and matrix (and interface or interphase), 
which by themselves are brittle or quasibrittle monolithic 
materials, the application of the Weibull distribution to describe 
damage initiation appears to be a more straightforward and 
reasonable exercise. Indeed this approach was taken by 
Guillaumat and Lamon (1996) and by Lamon et al. (1998) for 
their simulations of the matrix damage evolution of a woven 
CMC using their isotropic strength response unit-sphere model. 
For fiber (bundle) strength modeling in a PMC, Harlow and 
Phoenix (1978) indicate that WLT is an operative mechanism 
and that the strength distribution is approximately Weibull 
within the probability range of interest. This argument is also 
supported by Mahesh et al. (1999, 2002), Mahesh and Phoenix 
(2004), and Landis et al. (2000). The review article(s) by 
Nemeth and Bratton (2010, 2011) discuss in more detail some 
of the theoretical work and simulations that others have per- 
formed regarding strength scatter and size effect in composites. 

In this report, the unit-sphere methodology is described first 
for an isotropic material strength response followed by the 
modeling needed to enable the prediction of strength response 
for a transversely isotropic material in fast fracture (which is 
the intrinsic material strength without regard to time- or cycle- 
dependent strength degradation). These extensions are for 
(1) flaw-orientation anisotropy, whereby a preexisting micro- 
crack has a higher likelihood of being oriented in one direction 
over another direction, and (2) critical strength or fracture 
toughness anisotropy, whereby the level of critical strength a lr 
or fracture toughness K lc for mode I crack propagation 
changes with the orientation of the microstructure (for a 


tensile mode of failure only). These modeling extensions are 
then demonstrated for a (macroscopically) homogeneous 
material — in this case, nuclear-grade graphite — followed by 
an application to a heterogeneous material — a unidirectional 
glass/epoxy PMC. Only biaxial tension/compression loading 
was considered here. The intent is to show the progression of 
the model from a homogeneous isotropic material, to a homo- 
geneous anisotropic material, and then to a heterogeneous 
anisotropic composite material. This includes illustrating the 
consequence of shear sensitivity, global fracture planes, and 
failure modes on multiaxial strength. 

The unit-sphere methodology is an attempt to provide a 
mechanistic basis to the problem of predicting the strength 
response of a composite under multiaxial loading that is 
improved in comparison to polynomial interaction equation 
formulations such as Tsai-Wu, Tsai-Hill, and Hashin, among 
others. The unit-sphere methodology is not a true physics-based 
simulation whereby cracks nucleate or initiate from preexisting 
flaws, grow, and interact to progressively damage the material. 
However, those methods are computationally intensive and do 
not establish the statistics of failure without further considera- 
tion of the size distribution, the spatial distribution, and the 
orientation biases of inherent source flaws or defects. 

2.0 Model Description 

2.1 Unit-Sphere Model for Isotropic 
Response 

In the Batdorf unit-sphere theory, the incremental failure 
probability APjy under an applied multiaxial state of stress 2 at 
a given location in the component can be described as the 
product of two probabilities: 

AP /F (l,a b9C ,AL)=AP 1( ,P 2( , (1) 

where A P xv is the probability of the existence in the volume 
AL of a crack having an equivalent critical strength between 
e>\ e qc and a leqc + A a leqc . Critical strength o Ic is defined as the 
remote, uniaxial fracture strength of a given crack in pure 
mode I loading. The term <3\ eqc denotes an effective 
(or equivalent) critical mode I stress from applied multiaxial 
stresses. The second probability, P 2 v, denotes the probability 
that a crack with critical strength 0 \ eqc will be oriented in a 
direction such that an effective stress 0 \ eq (which is a function 
of fracture criterion, stress state, and crack configuration) 
satisfies the condition 0 \ eq > c>\,. qc . An incremental volume A V 
is used in Equation (1) because an infinitesimal volume dV 
cannot enclose a crack of critical strength c Ie?c and associated 
critical crack length of a c . 
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The effective stress 0 \ eq represents an equivalent normal 
stress on the crack face from the combined action of the 
normal stress g„ and the shear stress x. The microcrack 
orientation is defined by the angular coordinates, a and 
where the direction normal to the plane of the microcrack is 
specified by the radial line defined by a and [3 (see Fig. 1(b)). 
For the sake of brevity, the development of the effective stress 
equations is not shown (for details, see Nemeth et al. (2003, 
2005)). For a penny-shaped crack with the Shetty mixed-mode 
fracture criterion (Shetty (1987)), the effective stress becomes 


*3\eq 


C>„ + , / CT„ + 


4t 


|_C(2 - v)J 


(2) 


where v is Poisson’s ratio and C is the Shetty shear- 
sensitivity coefficient, with values typically in the range 
0.80 < C< 2.0. As C increases, the response becomes 
progressively more shear insensitive. Shear increases the 
equivalent stress as shown in Equation (2), and this has a 
deleterious effect on the predicted material strength. For a 
penny-shaped crack with a material having a Poisson’s ratio of 
about 0.22 and C =0.80, 0.85, 1.05, and 1.10; Equation (2) 
approximates, respectively, the following criteria: Ichikawa’s 
maximum energy-release-rate approximation (Ichikawa 
(1991)), the maximum tangential stress (Erdogan and Sih 
(1963)), Hellen and Blackburn’s (1975) maximum strain- 
energy-release-rate formulation, and colinear crack extension. 
The value of C also can be fit empirically to experimental 
data — either on introduced cracks (as was done in Shetty 
(1987)) or on specimens being tested multiaxially. 

The strength of a component containing a flaw population is 
related to the critical flaw size, which is implicitly used in 
Batdorf s theory. Batdorf and Crose (1974) describe APi^as 


AD . J. ^0 V ( CT I eqc ) 

A P w = AV — -da 


da 


I eqc 


I eqc 


and P 2 v is expressed as 


4n 


( 3 ) 


( 4 ) 


where i p{a Ie?c ) is the Batdorf crack-density function and 
11(2, Oi eqc ) is the area of the solid angle projected onto the unit 
radius sphere in stress space (see Fig. 1) containing all the crack 
orientations for which <j\ eq > 0 \ eqc for the applied far-field 
multiaxial stress state 2. The infinitesimal area, <14, on the unit 
sphere represents a particular flaw orientation (a direction 
normal to the flaw plane), and 0 \ eq is an equivalent, or effective, 
stress, which is a function of an assumed crack shape 


x 



«x 



Figure 1 . — Projection of equivalent stress onto a unit radius 
sphere, (a) Cauchy stress components on an infinitesimal 
tetrahedron resolving a normal stress, a„, and a resultant 
shear stress, x, on a plane normal to the direction defined by 
angular coordinates a and p. (b) Projection of equivalent 
stress onto a unit radius sphere in the global coordinate 
system. The unit radius sphere represents all possible flaw 
orientations, where ai eq is an equivalent (or an effective) 
stress, a x , a y , and a z are normal orthogonal stress compo- 
nents, and ixy, Tyz, and Tzx are shear stress components. 

An infinitesimal area, d/4, on the unit sphere represents a 
particular flaw orientation (a direction normal to the flaw 
plane), and ai eq is a function of an assumed crack shape and 
multiaxial fracture criterion. 
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and multiaxial fracture criterion. The constant 4k is the 
surface area of a unit radius sphere and corresponds to a solid 
angle containing all possible flaw orientations. 

The probability of failure for the arbitrarily loaded volume 

Vis 


Pfy = 1 - exp 



O^o^dr^CTi 

eqc ) 

4k dcr Ia/c 



(5) 


where a Ic9iinax is the maximum effective stress that a flaw 
could experience from the given stress state 2 and for all 
possible flaw orientations. The crack-density function i"l({ci|„ /r ) 
is independent of the stress state and has been approximated 
by a power function (Batdorf and Heinisch (1978)). This leads 
to the Batdorf crack-density function of the form 

I \ k BV oZ r (x, y,z,a,P) 

Wl e q c)= (6) 

GoV 

where x, y, and z correspond to the location, a and (1 are 
orientation angles, m v is the shape parameter, and kgy is the 
normalized Batdorf crack-density coefficient. 

The scale parameter a a y corresponds to the stress level 
where 63.21 percent of tensile specimens with unit volumes 
would fracture. A similar quantity, the characteristic strength 
Ogy, is the stress level (which is the peak stress for a nonuni- 
formly stressed specimen — such as a specimen undergoing 
four-point bending) at which 63.21 percent of specimens fail. 
The characteristic strength is functionally related to the scale 

parameter by cjgy = a oV jv^ mr as explained in other publica- 
tions, such as Nemeth et al. (2003, 2005). The effective 
volume V e is an equivalent amount of volume under uniform 
loading (with no stress gradients) that produces the same 
probability of failure as a specimen of volume V under loading 
where stress gradients may or may not be present. For a pure 
tensile specimen where stress distribution is uniform in the 
gauge section, the effective volume is equivalent to the 
volume of the gauge section. When stress gradients are 
present, a more complicated integral relationship is required to 
determine V e (Nemeth et al. (2003, 2005)). 

The scale parameter o oV has dimensions of stress x 
(volume ) ]/mv , where my is the shape parameter (Weibull modu- 
lus) — a dimensionless parameter that measures the degree of 
strength variability. As my increases, the dispersion decreases. 


The normalized Batdorf crack-density coefficient kgy is a 
function of the microcrack geometry and mixed-mode fracture 
criterion chosen. The value of kgy normalizes Equation (5) to 
the two-parameter Weibull equation for a uniaxial stress state 
(see Nemeth et al. (2003, 2005)). It also means that 
Equation (5) takes on the characteristics of a two-parameter 
Weibull distribution. The Weibull distribution is obtained if 
the defect size distribution is described by a power law. A 
Gumbel-type extreme-value distribution should be obtained if 
the defect size is exponentially distributed. Typically brittle 
material strength is described with the Weibull distribution. 
The normalized Batdorf crack-density coefficient kgy, the 
scale parameter a oV , and the Weibull modulus m v are 
evaluated from experimental inert strength (fast-fracture) data 
of specimens. 

To determine the probability of failure over V, one must 
evaluate P 2 y (in Eq. (5)) for each elemental volume AF, within 
which a uniform multiaxial stress state ct is assumed. The solid 
angle Q(2, a Ie?c ) depends on the selected fracture criterion, the 
crack configuration, and the applied stress state. For multiaxial 
stress states, with few exceptions, Q(2, <Ji eqc ) must be deter- 
mined numerically. Fortran was used to develop the unit- 
sphere numerical algorithms for this work. 

For a sphere of unit radius (see Fig. 1), an elemental surface 
area of the sphere is cL4 = sin a d(3 da. If we project onto the 
spherical surface the equivalent (effective) stress a k(? (2, a, (3), 
the solid angle £2(2, Oi eqc ) will be the area of the unit sphere 
containing all the projected equivalent stresses satisfying 
Gi eq > <3\ eqc . Because of the symmetry of the stresses, this 
integration can be performed over one -half of the sphere 
(0 < a < 7t/2); therefore, 

n ^4^ = (^) l 0 2 7o V2// ^’ 0l -) Sina d<X dP (7) 

where 


H\G\eq i ®Ieqc > 

1=1 

®leq — ®leqc 

H I eq 5 ^ I eqc ) 

1=0 

®Ieq ^ ®leqc 


and H is the Heaviside function. Substituting Equations (6) and 
(7) into Equation (5) and integrating with respect to <3\ eqc results 
in (see Batdorf (1978a, b) and Nemeth et al. (2003, 2005)): 


Pjv = 1 - exp 


k B y f |- 2 ^f V 2 fa Ie? (x,y,z,a,p 


2k 


U ZTt j- 71/ 
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For a given incremental volume, a i eq (x, y, z, a, P) is the 
projected equivalent stress over the unit radius sphere in stress 
space as depicted in Figure 1. Using the power-law crack- 
density function r\ v of Equation (6) will characterize the 
Batdorf model in Equation (8) as a form of the Weibull 
distribution. Equation (8) circumvents the more involved 
numerical integration of Q (Z, <5\ eqc ). 

It should be mentioned that, in the numerical implementation 
of Equation (8) when g„ becomes compressive, its value is set 
to zero in Equation (2) and t is adjusted by subtracting the 
absolute value of g„ from it. This happens when the absolute 
value of (g„/t) < 0.5 in the algorithm so that the adjusted value 
of t is never allowed to be less than zero. When the absolute 
value of (g„/t) > 0.5 (which was a somewhat arbitrarily set 
value from trial and error experience), then 0 \ eq is set to zero. 
These conditions are imposed to avoid an abrupt transition so 
that, when g„ goes from tension to compression, there will 
never be a failure of the flaw (regardless of the value of i). By 
subtracting the absolute value of a„ from i, the effective stress 
a Ieq is reduced. This provides a penalty (in a simplistic way) for 
the effect of friction on the crack face by increasing the 
resistance to mode II and III failure. This was found to be a 
practical way to avoid abrupt transitions in failure probability 
tendencies when going from tensile to compressive loading 
(such as seeing discontinuities in failure envelopes when 
contour lines of constant probability of failure are being 
plotted). 

2.2 Unit-Sphere Model Extended for Trans- 
versely Isotropic Strength Response 

Two different physical mechanisms were considered in 
order to extend the unit-sphere model to account for aniso- 
tropic strength response: (1) flaw-orientation anisotropy, 
whereby a preexisting microcrack has a higher likelihood of 
being oriented in one direction over another direction, and 
(2) critical strength or fracture toughness anisotropy, whereby 
the level of critical strength Gi c or fracture toughness K\ c for 
mode I crack propagation changes with regard to the orienta- 
tion of the microstructure (for a tensile mode of failure only in 
this case). Flaw-orientation anisotropy was previously consid- 
ered by Buch et al. (1977), and critical strength anisotropy was 
previously considered by Batdorf (1973). Both models were 
developed to model monolithic graphite anisotropic strength 
response. Extensions to these models are described herein for 
a transversely isotropic strength response. The extensions 
include shear sensitivity for flaws, an improved functional 
form for the anisotropy equations, and consideration of 
multiple failure modes. A corollary of the unit-sphere model is 


that the probability density distribution of the orientation of 
critical flaws in a multiaxial stress state can be obtained. This 
is described separately in Nemeth (2013). 

2.2.1 Flaw-Orientation Anisotropy 

Flaw-orientation anisotropy refers to the situation where a 
flaw has a higher likelihood of being oriented in one direction 
than another for a given critical strength. Thus, a material will 
be stronger on average in one direction than another. An 
isotropic brittle material is equally strong in any direction, and 
thus its flaws are uniformly randomly oriented. However, for 
components made by processes such as extrusion or hot- 
pressing, which induce texture, a bias will exist in the orienta- 
tion distribution of processing flaws. Also, components 
finished by surface grinding will contain machining damage in 
the form of surface cracks that are oriented parallel and 
transverse to the grinding direction. For composite materials, 
the interface, or an interfacial layer between the fiber and the 
matrix, can act as a flaw with an orientation bias that induces 
an anisotropic strength response. 

For volume-distributed flaws, flaw-orientation anisotropy 
relative to a material coordinate system is modeled by intro- 
ducing a probability density distribution p into the unit- 
sphere formulation (see Buch et al. (1977). Then for P :v in 
Equation (4), 

Pir= J J £>( a ,P )H(<j le q ,Gi £ , ?c )sina da dp (9) 

where 

p(a,P) = 2n n , 2 ^ a ^ ( 10 ) 

[ [ ~^(a,p) sina da dp 

J 0 J 0 

and 

P y^lsq •> ®leqc ) — ^ ^1 eq — ®]eqc 

H ( *5\eq ? ^ 1 eqc ) — ^ ®leq < ® leqc 

The term p(a, P) da dp denotes the probability that a crack 
of critical strength a ieqc is oriented in the range between a and 
(a + da) and between P and (P + dp), where the orientation of 
the microcrack is described by the vector that is normal to the 
plane of the microcrack. The function £^(a, P) describes the 
degree of anisotropy of flaw orientation, where the normal 
direction to the flaw plane is given by angles a and p. 
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Equations (9) and (10) modify Equation (8) for probability along the a x direction; and for the equatorial-belt distribution, 
of failure as it is evaluated in the o y (or optionally the a z ) direction. 


P fV =1 - exp 


- r r 2 * r it / 2 / / ( 

J, Jo Jo 


°i e<7 ( x ,y> z ,wP 


sin a da dp d V 


( 11 ) 


where the normalized crack-density coefficient k BV is 
evaluated numerically for a uniaxial stress state applied in one 
of the material coordinate system axis directions. 

For a transversely isotropic strength response, C, in 
Equation (10) is only a function of a. Buch et al. (1977) 
introduced a cosine power function for C(a) = [cos (a)]'*’ 
where <j) is a constant. This relation is modified herein to 
enhance the functional flexibility: 


?(“) = 

q(a) = 0 
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( 12 ) 


The separate polar-cap and equatorial-belt distributions are 
introduced to describe individual, and distinctly different, 
failure modes. For a unidirectional fiber-reinforced composite, 
the polar cap can be used to represent the fiber strength 
distribution and the equatorial belt can be used to represent the 
matrix-fiber interface. The equatorial-belt and the polar-cap 
distributions can be considered to equivalently represent 
global failure planes, which are referred to as “action planes” 
in the Puck multiaxial strength model for composites (see Lutz 
(2006) for a description). The angular width of the belt or cap 
with regard to angle a indicates the maximum extent of scatter 
of a flaw orientation. Flaws or action planes oriented outside 
of the scatter bands of the cap and belt are assumed not to 
exist or to not contribute to the likelihood of failure. 

2.2.2 Critical Strength or Stress Intensity Factor 
Anisotropy 


Alternatively, 

q(a) = 0 


f 

q(a) = 

V 

■P ^ 

tor — A t < a < — 

2 2 J 

where A and 4> are constants (with subscripts L or T) that 
control the degree of anisotropy, and 0 < A < ti/ 2 and (j) > 0. 
When A = n/2 and <|) = 0, an isotropic strength response is 
obtained. Equations (12) and (13) are defined for one-half of 
the unit sphere (the top half is as shown in Fig. 1, where 
0 < a < n/2). Referring to Figure 2, Equation (12) represents 
the “polar-cap,” or longitudinal, L distribution of flaws and 
Equation (13) represents an “equatorial-belt,” or transverse, T 
distribution of flaws. The polar-cap distribution describes 
crack planes symmetrically distributed (centered) about a 
plane (in this case the o y -o / plane), and the equatorial-belt 
distribution describes crack planes symmetrically distributed 
(centered) along a line (in this case, the cr x axis). For the 

polar-cap distribution, k BV is evaluated for a uniaxial stress 


for 0<a< A r 
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(13) 


Batdorf (1973) approached strength anisotropy using the c>\, 
strength ellipsoid approach (describing an ellipsoid rather than 
a unit sphere). In this report, strength anisotropy is used with 
the critical stress intensity factor (fracture toughness) K lc 
varying with the orientation angle on the unit sphere. This is 
functionally equivalent to o Ic varying with orientation, and it 
is also functionally equivalent to the size of the flaw changing 
with the orientation angle. For a CMC, where failure from 
loading in the fiber direction is by matrix cracking with large- 
scale fiber bridging, fracture toughness cannot be defined on 
the global scale of the structure. In that case, one has to 
alternatively consider critical strength as a metric. However, it 
is acceptable to use fracture toughness on the local scale at the 
crack tip, where micromechanics can account for the bridging 
explicitly. In this report, fiber bridging is not considered 
directly. Of first-order importance herein is that the local 
fracture toughness could change with the orientation of the 
flaw plane (or action plane) and the applied loading. The 
specific micromechanics of how this might occur is not 
considered here. 

The modeling approach taken here is similar to that 
described previously for flaw-orientation anisotropy. The criti- 
cal strength Oi, is defined as the fracture strength of the crack 
in mode I loading and is proportional to K lc . Therefore, for 

anisotropic K lc ( a, P) = (c Oi c (a, P)) = [ca ICjmax / If .(a,p)] , where 
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Figure 2. — Unit sphere with probability density distribution functions (PDFs) describing the anisotropy of the flaw 
orientation, where a and p are angular coordinates and A l and Arare constants in the flaw-orientation aniso- 
tropy function representing longitudinal and transverse distributions. Orientation is described with the normal 
to the crack plane. In this figure, two orientation functions are described: (1) a polar-cap distribution (Eq. (12)) 
describing crack planes symmetrically distributed (centered) about a plane (in this case the a y -a z stress 
plane) and (2) an equatorial-belt distribution (Eq. 13)) where crack planes are symmetrically distributed (cen- 
tered) along a line (in this case the a x axis). 


c is a constant (c = Y^fa^ with crack-shape geometry factor Y 
and critical crack length a c ), 0 | c max is the maximum value of 
<3\ c over the unit sphere (for all a and |3), and f\ c (a,P) is a 
normalized function expressing the degree of this anisotropy. 
The unit-sphere formulation described previously is used, 
except the Heaviside step function in Equations (7) and (10) is 
modified as follows 
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where <3\ eqc is substituted for Oi, to indicate generalized mixed- 
mode loading. 

For a transversely isotropic response where anisotropy is 
only a function of angle a, the function fj c { a) over the top 
half of the unit sphere (0 < a < 7i/2) is arbitrarily defined for 
the L distribution as 


/ic( a ) = 1 - 

for 0 < a < t, L 


f n i 
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and for the T distribution 


A. (a) = 1 for 
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(16) 


where ^ c T . y % It, >~l, and r T are constants. The L subscript 
relates to the polar-cap strength anisotropy distribution for 
crack planes symmetrically distributed (centered) about a 
plane (in this case, the cr y -C7 z plane), and the T subscript 
relates to the equatorial-belt strength anisotropy distribution 
for crack planes symmetrically distributed (centered) along a 
line (in this case the g x axis) (see also Fig. 2 for a reference 
frame). Figure 3 shows a schematic of Equations (15) and (16) 
in stress-strength space. 

Equations (12), (13), (15), and (16) are arbitrary functions 
chosen for their flexibility to fit to data. Other functions could 
be easily substituted. With this model, failure probability can 
be expressed in its most general form as 


a lc,x 

t 



Figure 3. — Schematic of Equations (15) and (16) in normalized 
strength space. Equation (15) correlates to the polar-cap 
distribution, and Equation (16) correlates to the equatorial- 
belt distribution, as shown in Figure 2. Flere, 'Ll, Z,t, yu and y r 
are constants in the critical mode I stress intensity anisotropy 
function representing longitudinal and transverse distribu- 
tions; f lc (a) is the normalized anisotropy function of the 
critical mode I stress-intensity factor, K\ c , or critical mode I 
strength, cn c , as a function of angle a; and ct, (; x , a, cy , and 
CT, (;Z are the normalized (by ai c , m ax) orthogonal critical 
strength components. 

example, the normalizing Equations (15) and (16) to modify 
m y (dividing m v by A (a) ). 

Experimental data of K lc varying with orientation mapped 
from indentation testing with the indenter oriented at various 
angles can hypothetically be used to describe the anisotropic 
strength response. Figure 4 shows an example from Corbin 
et al. (1988) using Equation (15) to fit to experimental data for 
anisotropic K ic measured from a Knoop indenter. Note that a 
more recent journal article by Quinn and Bradt (2007) recom- 
mended that indentation testing no longer be used for fracture 
toughness testing. However, in this case, only the relative 
difference of Aj, from one orientation to another is required 


P/y = 1 - exp 


- (• r 2 ti f s/2 , /a, (x,y,z,a,(3 


y»r 


sin a da dpdK 


(17) 


For an isotropic distribution of the orientation of the flaws, 
p(a, P) is 1/271, which is consistent with Equation (8). 
Although not shown herein, the Weibull modulus can also be 
made anisotropic, varying with angles a and P using, for 


rather than an accurate value for K Ic . The determined 
functional relationship in conjunction with the unit-sphere 
approach enables the stochastic strength response to be 
predicted from the multiaxial stresses with this model. 
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Figure 4. — Indentation testing used to determine the relation- 
ship between the critical mode I stress intensity factor, K| C , 
and the whisker axis direction. Here the constants in the K\ c 
anisotropy function representing the longitudinal distribu- 
tion are ^ = n/2 radians (90°), y l = 1 .0, and n = 1 .23 from 
Equation (15), and a is an angular coordinate. Experimental 
results from Corbin et al. (1988). 


Alternatively, specimens could be cut at various angles to a 
given direction in order to obtain strength properties. From that 
information, the necessary functional relationship of strength 
with angle would also be obtained for model calibration. 

2.3 Failure Modes and Overall Probability of 
Failure 

A material is assumed to have n individual failure modes. 
Material reliability (probability of survival P s where P s = 1 - Pj) 
is assumed to be the product of the survival probability of all of 
the failure modes: 


^=fl fe) 

i=l 

where each failure mode (denoted with the subscript i) has its 
own unique failure criterion and parameter description. For an 
isotropic material volume, two failure modes can be assumed: 
one for tensile failure and the other for compressive failure. 
The tensile-failure mode is described by Equation (8) and the 
effective stress relation of Equation (2). The compressive- 
failure criterion is assumed to be of a different nature than the 
tensile-failure mode — possibly involving the interaction of 
arrested cracks (flaws with initial crack growth that subse- 
quently stops because of the angle of the crack growth and the 
interaction with the local stress field). Regardless, the Weibull 
distribution is assumed to describe the stochastic strength 
response phenomenologically. This is argued further in 
Nemeth and Bratton (2011) on the basis of the statistical 
uncertainty in experimental data and a simple version of 
Freudenthal’s (1968) uniform defect model (Nemeth and 
Bratton (2011), App. A). The compressive-failure mode is 
assumed to be controlled by shear stress, and a simple 
Tresca-like effective stress relation can be prescribed: 


®leq = 2l (19) 

The multiplier of 2 in Equation (19) is chosen so that the 
maximum effective stress on the unit sphere in pure uniaxial 
compression is equal to the applied compressive stress. When 
the normal stress component o„ on the crack face is tensile, then 
the value of 0 \ eq in Equation (19) is set to zero. In this manner, 
the integration about the unit sphere in Equation (8) proceeds. 
This methodology can be used with the flaw-orientation 
anisotropy model of Equations (9) to (13). Flowever, a version 
analogous to the critical strength anisotropy model of 
Equations (14) to (17) has not been developed. Therefore, for 
this report, an isotropic strength response was assumed with 
this failure mode, unless it was specifically stated that 
flaw -orientation anisotropy was also used with this model. 

Modeling compressive failure with a Weibull distribution is 
a phenomenological assumption here. Scatter in strength can 
be substantial in compression as exemplified in Kittl and 
Aldunate (1983) for over 500 cement cylinders. They found 
that the normal, lognormal, and three-parameter Weibull 
distribution could not be chosen conclusively over each other 
as the best fit to the data. In a more rigorous treatment, Alpa 
(1984) extended the Batdorf unit-sphere methodology to 
account for compression by including the frictional effects of 
the opposed crack surfaces in contact for the multiaxial 
loading of an isotropic brittle material. Alpa’s treatment 
required that the Weibull modulus in tension and compression 
were the same. Here, failure in compression is treated as a 
separate failure mode, which allows the Weibull modulus for 
the compressive strength response to be different than the 
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Weibull modulus in tension. Puck and Shurmann (1998) 
(as well as Pinho et al. (2005)) consider frictional effects in 
their solution of the fracture plane. 

Multiaxial strength models are typically calibrated to exper- 
imental rupture data for simple loadings representing particu- 
lar failure modes, such as uniaxial tension and uniaxial 
compression. Specimen geometries such as cylinders, beams, 
or circular disks loaded in tension, flexure, or compression are 
usually used, and repeats of each test are required in order to 
obtain the fracture statistics (and hence, the Weibull param- 
eters). For an isotropic material, stochastic strength data in 
tension and compression can be obtained without regard to 
orientation relative to a material axis. 

2.3.1 Unidirectional Composite Failure Modes 

For an anisotropic material, stochastic strength data must be 
obtained relative to the different material axes. For a unidirec- 
tional PMC or CMC, this means strength testing in tension and 
compression for loadings parallel and transverse to the fiber 
direction as well as shear testing parallel and transverse to the 
fiber direction. Four failure modes are assumed herein for a 
unidirectional composite: (1) fiber fracture failure mode in 
longitudinal tension (loading in the fiber direction), (2) failure 
from longitudinal compression, (3) failure from transverse 
tension (loading perpendicular to the fiber direction), and 
(4) failure from transverse compression loading. These are 
assumed to be noninteracting competing failure modes 
possibly confined (or coalescing) along specific macroscopic 
fracture planes. This approach is consistent with Puck’s theory 
that the failure of a composite is controlled or confined to 
various action planes. 

Tensile and compressive rupture loadings are used to esti- 
mate Weibull parameters, whereas shear loadings help to 
establish C. Failure from pure shear loading is not considered 
to be a separate failure mode here because it is coupled to 
tensile failure by the assumption that the same flaws that 
control the tensile-failure response also control the shear 
strength failure response. This is indicated in Equation (2), 
where mixed-mode loadings are coupled through an effective 
stress relation. This also implies that the Weibull modulus in 
transverse tension and shear loading (both matrix-controlled 
failure modes) are the same. Flowever, in the PMC review by 
Wisnom (1999), the Weibull modulus seems to be larger in 
shear tests than for matrix transverse tension. It was not clear 
whether or not this was an artifact of the testing and the 
complex failure mode that results. Regardless, the assumption 
that shear mode failure has the same Weibull modulus of the 
matrix transverse tension failure mode is conservative. 

For transverse compression, the shear loading on the 
unfavorably oriented flaws in the matrix are assumed to drive 
the failure response (see Eq. (19)). This is allowed herein to be 
a phenomenologically different failure mode from the tensile- 
failure mode of Equation (2). The Weibull modulus can 


therefore be different from that of the matrix tensile-failure 
mode as discussed in Section 2.3. 

Failure from transverse tension and failure from transverse 
compression loading are assumed to be matrix-controlled 
failure modes. This is handled by the unit-sphere model as 
described in Section 2.0. The fiber fracture failure mode for 
the fiber bundle is also assumed to behave in an approximately 
Weibull manner as discussed in Section 1.0. Flere the same 
unit-sphere model will be used; however, fiber fracture will be 
assumed to be confined (or localized) to a specific fracture 
plane normal to the fiber (the polar-cap distribution as shown 
in Fig. 2). This essentially means that fibers are assumed to 
fail only by uniaxial tensile loading. This is not to say that a 
more general three-dimensional unit-sphere model could be 
used here to account for multiaxial loading and possible 
strength anisotropy in the fiber. Flowever, for the demonstra- 
tion problem that follows, this simplified failure criterion 
suffices. For PMCs, the review by Wisnom (1999) indicates 
that there is a clear size effect in the matrix for transverse 
tensile loading consistent with the Weibull theory. 

For the failure mode resulting from uniaxial compressive 
loading along the fiber axis (longitudinal loading), another 
simple criterion is assumed. In this case, failure of the compo- 
site occurs when the longitudinal compressive stress in the 
fiber exceeds a critical value. This critical value is assumed to 
be Weibull distributed. This criterion is provided here only so 
that closed failure envelopes can be shown. Failure from 
longitudinal compression is an area where research and model 
development is still evolving. It involves mechanisms such as 
microbuckling and fiber kinking. These are topics beyond the 
scope of this report. Pinho et al. (2005) provide a good 
overview and review of this topic. The review by Wisnom 
(1999) for PMCs does cite evidence for a size effect for this 
mode of loading, although the few studies performed made 
drawing firm conclusions difficult. 

3.0 Examples 

In the two examples that follow for nuclear-grade graphite 
and a unidirectional PMC, probability-of-failure strength 
envelopes were made for various combinations of biaxial 
loading. The nuclear-grade graphite material is treated as a 
continuum, whereas the PMC is a heterogeneous material 
system that requires the determination of the individual 
constituent micromechanical stress fields. The nuclear-grade 
graphite example is included to show the effect of shear 
sensitivity on flaws in comparison to its effect on unidirec- 
tional PMCs. 

3.1 Application to Nuclear-Grade Graphite 

Figure 5 presents the predicted biaxial failure envelopes for a 
near isotropic IG-110 (Sookdeo et al. (2008)) and transversely 
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isotropic H-451 (Burchell et al. (2007)) nuclear-grade graph- 
ites. Testing was performed on axially loaded and internally 
pressurized thick-walled hollow graphite cylinders. This 
imparted controlled axial and hoop stresses in the specimen 
gauge volumes. The manufacturing process produced a mild 
strength anisotropy in the IG-110 graphite billets and a stronger 
anisotropic response in the H-451 billets. The cylindrical 
specimens were cut from the billet stock so that axial specimen 


loading tested the stronger material direction and that hoop 
stresses from the internal pressurization tested the weaker 
material direction. Two separate failure modes for tension and 
compression were assumed, and the Weibull parameters used in 
the modeling were calibrated as such. Multiaxial strength 
response and the effect of the shear sensitivity parameter C 
(Eq. (2)) on the natural flaws were predicted. 



Figure 5. — Example of predicted biaxial failure envelopes for a near isotropic and 
transversely isotropic nuclear-grade graphite. Two separate failure modes for 
tension and compression were assumed. The envelope of the predicted effect of 
the shear sensitivity parameter C (Eq. (2)) on the natural flaws is shown (with 
C ranging between 0.82 and 100). (a) Isotropic strength response predicted for 
50-percent probability of failure, Pf, for IG-1 1 0 grade graphite (Sookdeo et al. 
(2008)). (b) Transversely isotropic strength response predicted for H-451 
graphite (Burchell et al. (2007)) for various failure probabilities and values of C . 
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Figure 5(a) shows the predicted 50-percent probability-of- 
failure response for IG-110 grade graphite, where it was 
assumed that the strength behavior was isotropic (the actual 
material behavior was near isotropic). For the tensile response, 
the employed parameters are m v = 22.43, Cq V = 26.38, and 
v = 0.22 (assumed), and C ranges from 0.82 to 100. Note that 
a 0 i' is not specifically needed here because the stress state in 
the gauge section can be considered to be uniform. For the 
compressive response, m e y = 22.43 and <j% e v = 81.32 MPa 
were assumed because there were no available compressive 
strength data (subscript g denotes the compressive failure 
mode). A strongly shear-sensitive response was imposed with 
C= 0.82, and a shear insensitive response was imposed with 
C= 100. The model correlated best to the experimental data 
for C « 1.2, indicating a moderate shear sensitivity on the 
flaws. It should be mentioned that values for the Batdorf 
crack-density coefficient k BV are generally not provided in 
this report. These values are numerically calculated within the 
software from the experimental data and do not require any 
specific treatment by the user. This parameter is described 
further in Nemeth et al. (2003, 2005). 

Figure 5(b) models a transversely isotropic strength 
response for FM151 graphite, showing C ranging between 
0.82 and 100 and various failure probabilities. Anisotropy of 
K\ c with flaw orientation was assumed. For the tensile 
response, the parameters are m v = 6.58 and agy= 17.05 MPa in 
the axial direction, with v = 0.22. In the hoop direction, a char- 
acteristic tensile strength response of cJeK = 11.01 MPa was 
modeled, where the parameters for the equatorial-belt 
distribution in Equation (16) were ^ T = 1-571 radians (90°), 
y T = L0, and (C = 0.82; r T = 2.323), (C = 1.4; r T = 1.915), 
and (C = 100; r T = 1.708). The values for factor r r were found 
using a simple interval halving search algorithm for a user- 
specified level of failure probability and multiaxial stress state. 
The values for the parameters = 1.571 radians (90°) and 
y T = 1.0 were arbitrarily chosen on the basis of an assumption 
that the variation of K k: with orientation was gradual — not 
abrupt. This approach is consistent with that of Figure 4 from 
the results of Corbin et al. (1988) for an extruded silicon- 
nitride/silicon-carbide whisker system. For the FI-451 
graphite, no information was available on the variation of K ic 
with orientation nor was it known if (and to what degree) 
anisotropic residual stresses were present in this batch of 
material. Flowever, the assumption that Aj, would not vary 
abruptly with orientation appears to be reasonable. For the 
compressive strength response, m e y = 12.29 and Oq £V = 
54.39 MPa. An isotropic compressive strength response was 
assumed because there were no available compressive strength 
rupture data in the hoop direction. Also, as previously 
mentioned, an anisotropic version of this failure criterion for 


critical strength has not been developed. From Figure 5(b) the 
model that correlated best to the experimental data was for 
C « 1 .4, again indicating a moderate shear sensitivity on the 
flaws. 

Figure 5(a) is shown mainly for illustrative purposes. 
Models of FI-451 graphite assuming flaw-orientation aniso- 
tropy (Equations (9) to (13)) or that an anisotropic residual 
compressive stress was present were not attempted here. 
Nemeth et al. (2012) showed that the trends observed in the 
fracture data of an extruded FI-451 graphite log (data from 
Price (1976)) were consistent with an anisotropic compressive 
residual stress component. Flowever, any physical evidence to 
further support that conjecture was not available. 

It can be seen in both Figure 5(a) and (b) that the shear 
sensitivity on the flaws primarily affects quadrants two and 
four (the tension-compression quadrants) for the tensile-failure 
model using Equation (2). The effect of C on the first 
(tension-tension) quadrant is small and would be difficult to 
detect from experimental data. The failure envelope in 
Figure 5(a) for compressive loadings is qualitatively compara- 
ble to those of Alpa (1984). Although not shown, there are 
transition regions present where one failure mode becomes 
more dominant over the other. 

It should be noted that, although it is true that the param- 
eters of the model were tuned to the graphite data, the aniso- 
tropy parameters themselves do have physical meanings and 
they could, hypothetically, be tested independently. Also the 
model could be tested against other multiaxial loading scenar- 
ios. The fact that the model worked well in this case indicates 
that the modeling approach is viable. 

3.2 Application to a Unidirectional Polymer 
Matrix Composite 

The biaxial strength response for a PMC was predicted from 
the strength response of the individual material constituents. In 
order to achieve this, a finite element analysis (FEA) of a 
square-packed fiber-in-matrix unit cell was used to obtain the 
microstress distributions within the fiber and surrounding 
matrix under biaxial loading. The effect of residual stresses 
from material processing was also considered in the finite 
element (FE) unit-cell model. 

From the FEA of the unit cell, the maximum (either positive 
or negative) stressed points were identified in each material 
constituent. Various failure mode scenarios were then applied 
to these individual points under biaxial loading to estimate the 
overall failure probability response of the composite. This 
approach is analogous to the direct micromechanics method 
(DMM) of Sankar (see Zhu and Sankar (1998) and Sankar and 
Karkkainen (2003)) except that unit-sphere methodology was 
applied as the failure criterion. The assumed failure modes 
were the (1) matrix in tension, (2) matrix in compression, 
(3) fiber bundle in tension, and (4) fiber bundle in 
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compression. Three variations of the anisotropic unit-sphere 
model, designated as Ml, M2, and M3, regarding flaw and 
failure plane anisotropy, isotropic matrix strength, and 
anisotropic matrix strength, respectively, were examined. Best 
results (in comparison to the experimental data) were obtained 
by the M3 model that considered the matrix as a mildly 
anisotropic material as opposed to considering matrix failure 
as being confined to specific fracture planes or to considering 
the matrix strength as being isotropic. The responses to biaxial 
loading of the individual failure modes and overall composite 
failure probability are also shown. This example demonstrates 
how unit-sphere methodology can be used as a failure criterion 
within finite-element- or micromechanics-based software 
codes for brittle composite material constituents to predict the 
overall stochastic strength response. 

The unit-sphere methodology is intended for use with brittle 
or quasibrittle material constituents. For this example, rupture 
data for a unidirectional PMC material were used. The data 
come from Al-Khalil et al. (1996) and were part of the WWFE 
(Flinton et al. (2004)). These data have been used widely as a 
benchmark with which to compare new failure criteria. 
Because of that and because of the simplicity of the material 
system, specimen geometry, and loading, this particular 
WWFE exercise case was chosen. 

Note that a PMC behaves less brittle in shear than in uni- 
axial tension, showing evidence of damage accumulation 
through matrix microcracking prior to failure. Flowever, it is 
assumed here that the matrix material itself fails from shear in 
a brittlelike manner similar to brittle failure in tension. Flence, 
as currently conceived, the unit-sphere prediction for the shear 
failure mode in a PMC is indicative of the probability of 
damage initiation and not specifically of damage progression 
and ultimate failure. 

The PMC studied here was a unidirectional lamina under 
combined longitudinal and transverse loading. The material 
was an E-glass/epoxy composite with a fiber volume fraction 
of 60 percent. The rupture data are reported in Al-Khalil et al. 
(1996) and summarized and listed in Soden et al. (2004b). As 
reported in Soden et al., most of the results were obtained 
from the testing of nearly circumferentially filament-wound 
tubes under combined internal pressure and axial load. For the 
unit-sphere-based analysis that follows, the fibers were 
assumed to be oriented at 0° relative to the hoop direction; 
however, the actual fiber winding angle was ±5° from the 
hoop direction rather than 0°. The specimens had an inner 
diameter of 100 mm, were 300 mm long, and were between 
0.95 and 1.2 mm thick. The E-glass fiber reinforcement was 
Silenka 05 1L, 1200 tex; and the epoxy resin system was a 
Ciba-Geigy MY750/F1Y917/DY063 mixture. Soden et al. 
(2004a) provide the properties for the fiber, matrix, and the 
composite. The experimental failure stresses are the stresses at 
which leakage and fracture occurred. The stress-strain behav- 
ior was essentially linear elastic up to the point of failure for a 


circumferentially (pressure load only) stressed tube. In the 
biaxial tests, all fractures were within the gauge section, and 
these failures occurred without prior warning. 

This exercise examined the probability of failure response 
under biaxial loading of the fiber and matrix constituents at four 
highly stressed locations in an FE fiber-in-matrix micromechan- 
ics unit-cell model of the composite. The relative contributions 
of the individual failure modes, as well as the effects of stress 
concentration and the multiaxial stress state at each point 
location, were tracked. In addition, the effect of confining 
failure to “action planes” was examined. This exercise not only 
tested the capability of the unit-sphere model to predict the 
failure response of the composite, but it provided additional 
insight into the material’s mechanical failure behavior. The 
approach taken here of using a unit-cell model to determine the 
microstress distributions between the fiber and the matrix and to 
subsequently apply multiaxial failure criteria is consistent with 
the DMM of Zhu and Sankar (1998) and Sankar and Kark- 
kainen (2003) as mentioned in this section. 

Individual failure modes were assumed as described in 
Sections 2.3 and 2.3.1. The reliability of the composite was 
computed using Equation (18). It was assumed that component 
failure was confined to the gauge section of the composite 
tube specimen and that stresses were uniform throughout the 
gauge section (that is, the stresses did not vary along the 
(axial) length, circumference, and thickness of the gauge 
section for a given biaxial loading). Therefore, the effect of 
the test component geometry and loading on the stress 
distribution was not considered in the analysis (consistent with 
the WWFE). The stress-volume integration (the integration of 
the stress over the whole volume of the particular material 
constituent of the unit-cell model) of Equation (5) was not 
performed. Instead this analysis only considered the effect of 
worst-case stresses on composite failure probability, akin to a 
maximum-stressed-point failure criterion. Thus, the potential 
effect on failure probability of the effective volume changing 
with multiaxial stress state was not investigated. Strength for a 
given probability of failure will scale with V e Vmr , which 
means that the effect of changing V e on strength is of a lower 
level of sensitivity in comparison to a directly proportional 
response, and this sensitivity decreases as the value of m v 
increases. However, the benefit of this approach (of using a 
maximum-stressed-point failure criterion) is that only the 
effect of the stress state on the failure probability is considered 
for different multiaxial loading conditions, without the 
additional confounding effect of V e being included. This 
means that the comparison of the unit-sphere prediction to 
other failure criteria is more direct — only predictions based on 
stress state are compared. A goal of this report is to investigate 
the relative contribution of the highly stressed regions of the 
unit cell, considering location and failure mode, to the overall 
composite failure response. The effect of stress-volume 
integration of the unit cell will be considered for future work. 
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An FE model of a 3-by-3 square-packed arrangement of 
fibers embedded in a matrix was constructed to determine the 
peak stresses in the respective material constituents of the 
composite. Only the stresses in the center fiber-in-matrix unit 
cell were considered in the analysis. The larger 3-by-3 model 
was constructed so that boundary (edge) effects on the center 
unit-cell of the FE model would be minimized. Figure 6(a) 
shows the mesh for the 3-by-3 model, and Figure 6(b) shows 
the mesh for the center unit cell (at the midlength along the fiber 
axis) with an arbitrary tensile transverse and longitudinal 
multiaxial loading superimposed on the model (red is tensile 
and blue is compressive). The model was constructed in Abaqus 
(version 6.10 Dassault Systemes (2010)) with 135,540 C3D8 
linear brick elements with a 60-percent volume fraction of 
fibers. Symmetric boundary conditions were used that restricted 
out-of-plane motion for the left, bottom, and back faces in 
Figure 6(a). Kinematic coupling constraints from a point to a 
surface were applied about the opposite faces (right, top, and 
front) to constrain out-of-plane degrees of freedom (out-of- 
plane motion of a face was constrained to the motion of a single 
point). Transverse-to-the-fiber loading was applied using 
pressure loads applied to the right face (in the y direction) and 
axial-to-the-fiber loads were applied using pressure loads to the 
front face (in the x direction). Residual stresses resulting from 
processing were accounted for by imposing an initial 120 °C 
uniform temperature on the model and subsequently cooling to 
20 °C (120 °C is listed as the “stress free temperature” in Soden 
et al. (2004a)). Isotropic material properties were assumed for 
the fiber and the matrix, and these are listed in Table I (see 
Soden et al. (2004a) for a complete listing). Linear-elastic 
constitutive behavior was assumed. 

Soden et al. (2004b) and Al-Khalil et al. (1996) provide 
experimental rupture data of the composite. The strength values 
from uniaxial loadings, which were used to normalize the model 
to a 50-percent probability of failure, are shown in Table II. 
They are also listed in Soden et al. (2004b) as mean strengths. 
The biaxial loadings applied on the FE unit-cell model were in 
the range of the strength values shown in Table II. Various 
longitudinal (axial) loads were applied from -800 to 1400 MPa 
for fixed transverse loads of 40, 0, or -145 MPa. For these load 
combinations, the highest (in absolute value) stresses occurred 
in the regions indicated by the six points labeled in Figure 6(b). 
These points are designated henceforth as the “matrix left 
edge,” “matrix left interface,” “fiber left interface,” “matrix top 
edge,” “matrix top interface,” and “fiber top interface.” These 
peak stresses involved the o x , o y , and a z stress components, 
whereas the x xy , Xyz, and x zx shear stress components were 
generally negligible at these four points. 

Figure 7 shows the stress response in the matrix for a fixed 
applied transverse load of -100 MPa and a varying applied 
axial load between 0 and 1200 MPa without any residual 
stresses from processing (Figs. 7(a), (c), and (e)) and with 



(a) 



Left edge 
( 1 ) 

Left 

interface 
(2) and (3) 


Top edge (4) 


Top interface 
(5) and (6) 


Figure 6. — Finite element model of a portion of a unidirectional 
composite consisting of nine fibers in a square-packed 
arrangement, (a) Finite element mesh of 3-by-3 square- 
packed array of fibers embedded in a matrix, (b) Detailed 
mesh for center cell. A representative tensile transverse 
stress (o y applied left to right) and longitudinal stress 
(a x applied along the fiber axis) biaxial loading scenario is 
shown where red is tensile (middle of left and right edges) 
and blue is compressive (middle of top and bottom edges). 
The points shown with x’s correspond to highest stressed 
regions found for the biaxial loading cases examined. The 
left interface and top interface both include the matrix side 
and fiber side of the interface boundary. 
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TABLE I.— PROPERTIES USED FOR THE FIBER AND MATRIX COMPOSITE CONSTITUENTS 


Material 

Type 

Young’s modulus. 

Poisson’s 

Thennal expansion coefficient, 

constituent 


E, 

ratio, 

a „ 



GPa 

v 

lO^/C 0 

Fiber 

Silenka E-Glass 1200 tex 

74.00 

0.20 

4.9 

Matrix 

MY750/HY 9 1 7/DY 063 

3.35 

0.35 

58.0 


TABLE II.— COMPOSITE PROPERTIES 


Fiber 

volume 

fraction 

Longitudinal 
tensile strength, 
MPa 

Longitudinal compres- 
sive strength, 

MPa 

Transverse 
tensile strength, 
MPa 

Transverse compres- 
sive strength, 
MPa 

Stress-free 

temperature, 

C° 

0.6 

1280 a 

800 a 

40 a 

145' 

120 


^Assumed to be a median strength value. 


TABLE III.— RESIDUAL THERMAL STRESSES FROM 
PROCESSING AT THE VARIOUS POINT LOCATIONS 


Point location 

Axial stress, 
G x , 

MPa 

Transverse stress: 
y direction, 

Gy, 

MPa 

Transverse stress: 
z direction, 

G z , 

MPa 

Matrix left edge 

13.4 

-30.5 

20.2 

Matrix left interface 

12.4 

-28.0 

13.3 

Matrix top edge 

13.6 

20.1 

-30.4 

Matrix top interface 

12.3 

13.1 

-29.0 

Fiber left interface 

-18.9 

-28.1 

1.38 

Fiber top interface 

-18.9 

1.43 

-28.1 


superimposed residual stresses from processing (Figs. 7(b), 
(d), and (f)). This loading is over the region of the fourth 
quadrant for biaxial loading (applied transverse compression 
with longitudinal axial loading). All stress trends versus axial 
loading are linear — as would be expected for linear elastic 
analysis. Generally, the lowest stresses are for the transverse 
compressive load component a y for the left edge and left 
interface and the highest stresses are for axial load ct x at the 
top edge and top interface. For all stress components in 
Figure 7, the stresses between the top edge and top interface 
have similar values. This is also true for the left edge and left 
interface except in Figure 7(e) and (f) for the ct z stress compo- 
nent, where some divergence can be observed. Stress inter- 
actions are present but are relatively mild; nonetheless, they 
do affect the probability of failure response. For example, in 
Figure 7(c) for the left edge and left interface, increasing the 
axial load makes the a y stress more compressive. Figure 7 also 
shows that the thermal residual stresses (listed in Table III) 
have only a minor effect on overall stress magnitudes. 

In order to generate the multiaxial failure envelopes, a small 
Fortran program was created that extrapolated the <r x , o y , and o z 
stress components for six evaluation points (four for the matrix 
and two for the fiber, as shown in Fig. 6) for any combination of 
applied axial and transverse loads. This was done to avoid 
unnecessary repetition running the FE unit-cell model. Shear 
stress components were assumed to be negligible. Linear 
relationships such as those shown in Figure 7 were developed 
for constant applied transverse stresses of 40 and -100 MPa 
versus applied axial load (without thermal residual stresses). 
Then, for a given ratio of applied transverse stress to axial 


stress, the g x , g v , and ct z stress components for the six points 
(for 40-MPa transverse load in stress quadrants 1 and 2 or for 
-100-MPa transverse load in stress quadrants 3 and 4) were 
scaled to correspond to a 100-MPa load on the hypotenuse 
formed from the right triangle of the transverse and axial 
applied loads, which defined an angle 9 measured from the 
longitudinal loading axis. For various angles of 9, the ct x , a y , 
and cr z stress components could be calculated for the 100-MPa 
load defined from the hypotenuse. In this manner, the g x , c y , g z 
stresses were calculated versus angle 9, defining a circle with 
radius of 100 MPa in the applied biaxial load stress space. Spot 
checking with FEA was used to verify the Fortran program. 

So that the failure envelope could be generated, the stress 
levels along a constant angle 9 were proportionally scaled 
until a 50-percent probability of failure was obtained using the 
unit-sphere numerical algorithms, which were also written in 
Fortran. The scaling factor needed to achieve 50-percent 
probability of failure multiplied by the 100-MPa radius of the 
base circle determined the radial length from the origin for the 
given angle 9. Repeating this procedure for various angles of 9 
enabled the failure envelope to be constructed. The effects of 
residual thermal stresses were accounted for by simply adding 
the respective residual thermal stress ct x , ct v , and ct z compo- 
nents to the overall g x , o y , and o z stresses (which were 
previously multiplied by the scaling factor). The residual 
thermal stresses are a constant and are present regardless of 
the level of the applied biaxial loadings. 

The probability of failure calculation considered the four 
sampled points in the matrix and the two points sampled in the 
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Stress, a z , MPa Stress, o y , MPa Stress, o x , MPa 








Figure 7. — Stresses at four locations (shown in Fig. 6) in the matrix for a fixed applied transverse load of -100 MPa 
and a varying applied axial load, (a) a x stresses without thermal load, (b) a x stresses with superposed thermal 
load from processing, (c) a y stresses without thermal load, (d) a y stresses with superposed thermal load from 
processing, (e) a z stresses without thermal load, (f) a z stresses with superposed thermal load from processing. 
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fiber for all the failure modes. Equation (18) was used to 
determine the overall survival probability from each individual 
failure mode. For example, the highest effective stress on the 
unit-sphere for a given angle a and P of the four sampled points 
was determined for the matrix material. Integration over the 
surface of the unit sphere for all angles a and P provided a 
probability of survival for the four sampled points. This calcula- 
tion was performed independent of the volume associated with 
each sampled point. Applied in this manner, the methodology 
becomes a maximum-stressed-point failure criterion that also 
includes the effect of the orientation of the flaw. 

3.2.1 Biaxial Failure Envelopes for Three Assumed 
Models 

Various combinations of unit-sphere modeling assumptions 
were tried to examine their behavior relative to experimental 
data and to well-known phenomenological failure criteria such 
as Tsai-Wu and Tsai-Hill. The following subsections describe 
the modeling combinations. In all cases, a Weibull modulus of 
m v = 10.0 was arbitrarily assumed for all failure modes except 
as noted. Failure envelopes are shown for 50-percent probabil- 
ity of failure. A Weibull modulus of m v = 10.0 is reasonably 
representative of matrix failure modes. The fiber failure mode 
from longitudinal tensile loading usually has a higher Weibull 
modulus (this is corroborated by Wisnom (1999)). For this 
WWFE exercise from the Al-Khalil et al. (1996) experiments, 
the Weibull modulus for the longitudinally loaded specimens 
was estimated to be in the range of m v = 20.0 (from a sample 
of eight specimens). Flowever, the Weibull modulus from 
the matrix failure modes was unknown. Regardless, because 
the 50-percent probability of failure results are calibrated to 
the rupture strengths in Table II, the particular value of the 
Weibull modulus chosen had little visible effect on the 
presented results. It should be mentioned that Al-Khalil (1996) 
adjusted the reported rupture stresses to account for “addi- 
tional” transverse tensile stresses. This added some uncertain- 
ty as to the true nature of the failure modes of the 
longitudinally loaded specimens. 

Three model variations are described — Ml, M2, and M3 — 
for flaw and failure plane anisotropy, isotropic matrix 
strength, and anisotropic matrix strength, respectively. In 
model Ml, fracture is assumed to be confined to specific 
failure planes using the flaw-orientation anisotropy model of 
Section 2.2.1. In model M2, the failure response of the matrix 
is assumed to be isotropic (Section 2.1) and the failure 
response was not confined to any specific action planes. 
Model M3 is identical to model M2 except that the matrix 
strength response is assumed to be transversely isotropic using 
the model for critical strength anisotropy described in 
Section 2.2.2. These models and their subsequent results are 
described in the following subsections. 

3.2.2 Model Ml: Flaw/Failure Plane Anisotropy 

In model Ml (interactive and noninteractive versions), frac- 
ture is assumed to be confined to specific failure planes or 


“action planes” in a manner similar to Puck’s theory. The 
interactive version uses stress results from the unit-cell FEA, 
and the noninteractive version uses only the applied composite 
(macroscopic) stresses. The Ml model variation uses the flaw- 
orientation anisotropy model described by Equations (9) to (13). 

For the fiber failure mode in tension, the failure is assumed 
to result from microscopic flaws in the fiber volume that are 
highly aligned perpendicular to the fiber axis or are confined 
to global failure planes as mentioned in Section 2.3.1. Fiber 
failure, therefore, only becomes a function of the tensile stress 
applied along the fiber axis (assuming in this case that the 
shear traction on these planes is negligible). This corresponds 
to the polar-cap distribution as shown in Figure 2 (and 
Eq. (12)), whereby the fiber axis is parallel to the o x material 
coordinate system stress axis. Effective stress is calculated for 
this failure mode, using Equation (2). For the fiber compres- 
sion failure mode, only the compressive normal stress compo- 
nent is used to evaluate failure in the manner previously 
described in Section 2.3.1. Here A L is assumed to be 1.0° with 
a uniform distribution (equal likelihood) of orientation over 
that increment (4>^ = 0.0) for the polar-cap distribution. 

For the matrix failure modes of model Ml, the flaws or 
failure planes are assumed to be highly confined in an 
equatorial-belt distribution as depicted schematically in Figure 8 
(also Fig. 2) and described by Equation (13), whereby A r is 1.0° 
with a uniform distribution (equal likelihood) of orientation 
over that increment (4*7’ = 0.0). The tightly confined equatorial 
belt distribution is meant to approximate an interfacial flaw 
population (or global failure planes) between the fiber and the 
matrix. For the tensile-failure mode, Equation (2) is used to 
describe the effective stress, whereas for the compressive- 
failure mode, Equation (19) is used. 

Table IV lists the parameters for these various failure modes 
for the interactive Ml model without the effect of the residual 
thermal stresses from processing, and Table V lists the param- 
eters that include the effects of residual stresses from thermal 
processing. The characteristic strength <Jq V is set to the value 
where the probability of failure is 50 percent for the overall 
composite loads indicated in Table II. The values of the charac- 
teristic strength are a function of the stresses at the sampled 
points in the fiber and matrix for the unit-cell micromechanics 
FE model. These values reflect the levels of stress concentration 
in the fiber and matrix predicted from the FE model. 

A noninteractive version of the Ml model was also devel- 
oped. The noninteractive model has the same unit-sphere 
modeling assumptions as the interactive version of Ml except 
that it uses the global applied composite stresses — not the 
unit-cell microstresses of the fiber and matrix. The approach 
of using global applied composite stresses is similar to other 
polynomial failure criteria such as Tsai-Wu and Tsai-Hill. The 
noninteractive model did not consider the effect of thermal 
residual stresses. The parameters for the noninteractive Ml 
model are given in Table VI. The values for were deter- 
mined such that the composite strength values listed in 
Table II correspond to a 50-percent probability of failure. 
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Figure 8. — Model Ml — matrix fracture planes are confined to be randomly oriented normal to the fiber axis; 
a x , ay, and <j z , orthogonal stresses. 


TABLE IV.— INTERACTIVE MODEL Ml UNIT-SPHERE PARAMETERS WITHOUT 
RESIDUAL THERMAL STRESSES FROM PROCESSING 


[One-half angle of anisotro 

jy distribution, A, 1.0°; exponent of sine or cosine function, (|), 0.0.] 

Failure mode 

Model and mode a 

Weibull 

Characteristic 

Shetty 



modulus, 

strength, 

shear-sensitivity 



m y 

<5ev, 

constant, 




MPa 

C 

Fiber: tensile 

L flaw anisotropic 

10.0 

2147 

1.4 

Fiber: compressive 

L flaw anisotropic 

10.0 

1346 

NA b 

Matrix: transverse tensile 

T flaw anisotropic 

10.0 

79.19 

1.4 

Matrix: transverse compressive 

T flaw anisotropic 

10.0 

181.8 

NA b 


a Z, longitudinal; T, transverse. 
hNA, not applicable. 


TABLE V.— INTERACTIVE MODEL Ml UNIT-SPHERE PARAMETERS WITH 
RESIDUAL THERMAL STRESSES FROM PROCESSING 


[One-half angle of anisotropy distribution. A, 1.0°; exponent of sine or cosine function, (|>, 0.0.] 


Failure mode 

Model and mode a 

Weibull 

modulus, 

111 y 

Characteristic 

strength, 

C>9F, 

MPa 

Shetty 

shear-sensitivity 

constant, 

C 

Fiber: tensile 

L flaw anisotropic 

10.0 

2128 

1.4 

Fiber: compressive 

L flaw anisotropic 

10.0 

1358 

NA b 

Matrix: transverse tensile 

T flaw anisotropic 

10.0 

57.15 

1.4 

Matrix: transverse compressive 

T flaw anisotropic 

10.0 

232.9 

NA b 


a Z, longitudinal; T, transverse. 
b NA, not applicable. 
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TABLE VI.— NONINTERACTIVE MODEL Ml UNIT-SPHERE PARAMETERS 
[One-half angle of anisotropy distribution, A, 1.0°; exponent of sine or cosine function, (|>, 0.0.] 


Failure mode 

Model and mode a 

Weibull 

modulus, 

m v 

Characteristic 

strength, 

<J0F, 

MPa 

Shetty 

shear-sensitivity 

constant, 

C 

Fiber: tensile 

L flaw anisotropic 

10.0 

1328 

1.4 


L flaw anisotropic 

20.0 

1304 

1.4 

Fiber: compressive 

L flaw anisotropic 

10.0 

829.9 

NA 5 


L flaw anisotropic 

20.0 

814.8 

NA b 

Matrix: transverse tensile 

T flaw anisotropic 

10.0 

41.49 

1.4 


T flaw anisotropic 

20.0 

40.74 

1.4 

Matrix: transverse compressive 

T flaw anisotropic 

10.0 

150.4 

NA 5 


T flaw anisotropic 

20.0 

149.7 

NA b 


“L, longitudinal; T, transverse. 
'’NA, not applicable. 



O Experimental data 

Tsai-Wu 

Tsai-Hill 

Ml — Noninteractive; m\/= 10 

Ml — Noninteractive; m\/ = 20 

♦ Ml — Without thermal 

• Ml — With thermal 


Figure 9. — Comparison of failure models for 50-percent probability of failure for the combined 
failure modes of the Ml (anisotropic flaw plane/global failure plane) model relative to the 
macroscopic applied stresses. Also shown is a non interactive model and Tsai-Wu and Tsai-Hill 
models for comparison. The failure envelopes for the two interactive models for the Weibull 
modulus for the volume-flaw failure mode, mv = 10, are shown with and without considering the 
predicted residual stresses resulting from processing. For the matrix tensile failure mode, Shetty 
shear-sensitivity constant C =1.4. 


3.2.2. 1 Model Ml Results 

Biaxial failure envelopes for a predicted overall 50-percent 
probability of failure for the combined failure modes (via 
Eq. (18) of Section 2.3) were generated with the unit-sphere 
algorithm for the Ml models for the four failure modes. These 
results are shown in Figure 9 relative to the macroscopic 
(composite) applied stresses. The Shetty shear-sensitivity 
parameter C= 1.4 for the tensile-failure modes was used. This 
parameter only has a limited effect, which will be shown later 
for models M2 and M3. Results from the Tsai-Wu and Tsai-Hill 
failure criteria are also shown in Figure 9 for comparison. The 
Tsai-Wu and Tsai-Hill criteria, as used here, are based on global 


applied stresses and, therefore, do not consider the effect of 
thermal residual stresses. Figure 9 also shows the overall failure 
envelopes considering all failure modes for the noninteractive 
Ml model for m v = 10 and 20 (where all failure modes have the 
same Weibull modulus value for each failure envelope) compar- 
ing the effect of m v on the shape of the failure envelope. 

Figure 9 shows the adjusted experimental data (Al-Khalil 
et al. (1996)) used in the WWFE exercise. What is immedi- 
ately apparent is that the Tsai-Wu criterion fits to the (limited 
amount of) experimental data quite well. Tsai-Hill also fits the 
data well in the fourth quadrant (longitudinal tension and 
transverse compression) but is conservative in the first 
quadrant (axial tension and transverse tension). The Ml 
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models fit less well to the data in the fourth quadrant. The Ml 
models are significantly different from Tsai-Wu and Tsai-Hill 
in all quadrants. In fact the noninteractive Ml model is nearly 
identical to the maximum principal stress failure criterion (see 
Fig. 3 in Zhu and Sankar (1998), for example) except that the 
corners of the envelope are slightly rounded because of the 
effect of Equation (18) for separate failure modes and the 
value of Weibull modulus m v that was used. The effect of m v 
is shown for the noninteractive model where the higher value 
of my = 20 appears to sharpen the “corners” of the failure 
envelope in comparison to the envelope for m v = 10. The 
noninteractive Ml results are also qualitatively similar in 
appearance to those presented in Puck and Schumann (2004). 

The interactive Ml model does show some differences from 
the noninteractive model in Figure 9. This difference is most 
pronounced in the fourth quadrant, and it illustrates the effect 
that the microstresses have on the predicted failure probability. 
For the most part, the envelopes with and without the account- 
ing of thermal residual stresses from processing were quite 
similar, except in the first quadrant, where the results without 
the residual thermal stresses were noticeably unconservative. 
The reason for the overall similarity was that the failure envel- 
opes were normalized to the experimental test data of Table II 
(which already include any effect from thermal residual 
stresses). Interestingly, the failure envelope without thermal 
load bears a qualitative resemblance to results from Zhu and 
Sankar (1998) for the particular DMM case (Von-Mises failure 
criterion for the fiber and matrix and no interfacial failure 
criterion) that was illustrated in Figure 3 of that paper. Some of 
the reasons that the interactive Ml failure envelopes took the 
shapes that they did are examined next. This line of inquiry was 
not pursued in the Zhu and Sankar (1998) paper. 

3.2.2.2 For Ml: Tracking the Response of Individual 
Sampling Points and Failure Modes 

Figure 10 shows the interactive Ml model for the individual 
sample points at 50-percent probability of failure relative to the 
applied macroscopic composite stresses for both the case with 
and the case without the thermal residual stresses. The failure 
contour lines closest to the origin dominate the failure response 
and create the enclosed region about the origin that results in the 
overall failure envelope shown in Figure 9. What is immediately 
apparent in Figure 10 is that the left-edge sampling point 
dominates the matrix failure response followed by the left 
interface sampling point for matrix tensile failure. The matrix 
top edge and top interface sampling points had little influence. 

For the fiber failure modes (tension and compression) — the 
fiber left edge and fiber top edge — sampling points gave 
similar results as indicated by the vertical lines in Figure 10. 
These fiber failure modes indicate little influence from the 
matrix and transverse loading on fiber failure, as would be 
expected based on failure mode modeling assumptions and the 
small relative influence that the matrix has on the fibers. For 
the matrix, the influence of changing stress state is seen as a 


deviation from a horizontal line. For example, in the fourth 
quadrant there is an increasing compressive stress for a y for 
the left edge and interface sampling points with increasing 
longitudinal (axial) load as shown in Figures 7(c) and (d), 
which influences the matrix compressive-failure mode by 
increasing the probability of failure with the increasing 
longitudinal load. The decreasing compressive stress for c z 
with increasing longitudinal load shown in Figures 7(e) and (f) 
has little influence because the magnitude of the compressive 
stress a z is significantly smaller than that of cr y regardless of 
the level of applied longitudinal stress (within the range of 
interest). The o x stress had no effect on the matrix probability 
of failure (see Figs. 7(a) and (b)) because the Ml model 
assumed that the flaw planes or global failure planes were 
essentially parallel to that direction of stress. Finally, the 
differences between the individual-sampled-point probability- 
of-failure contour lines without thermal residual stresses 
(Fig. 10(a)) and with thermal residual stresses (Fig. 10(b)) are 
not that significant, except perhaps in the first quadrant. 
Figure 11 shows the combined effect of all the sample points 
on the individual failure modes for a 50-percent failure 
probability for both cases with and without the residual 
thermal stresses present. Figure 11 provides a clearer picture 
of the relative contribution of the individual failure modes. 

3.2.3 Model M2: Isotropic Matrix Material 

In model M2, the matrix was assumed to be an isotropic 
strength response material, and the compressive-fiber failure 
mode was not used. This model only used the unit-cell 
microstresses. Because the strength response of the matrix was 
assumed to be isotropic, the M2 model also predicts the matrix 
strength response for applied uniaxial tensile and compressive 
loading in the longitudinal (fiber) direction. This differs from 
Ml, where flaw planes or failure planes were oriented parallel 
to the fiber direction, and hence, matrix failure could not occur 
regardless of the magnitude of load (infinite matrix strength in 
this direction, meaning that the fiber was always assumed to 
fail before the matrix). 

M2 used the same tensile-stress fiber-failure-mode model 
as Ml. The matrix tensile-failure mode and the matrix 
compressive-failure mode both used an isotropic strength 
response model of the unit sphere (Section 2.1). The effective 
stress equations (Eqs. (2) and (19)) were used for the tensile and 
compressive matrix failure modes, respectively. Table Vll lists 
the parameters for the failure modes without the effect of the 
residual thermal stresses from processing, and Table VIII lists 
the parameters that include the effects of residual stresses from 
thermal processing. Only results for the Shetty shear-sensitive 
parameter C = 1.4 for the tensile-failure modes are given. As 
described in Section 3.2.2 for model Ml, the characteristic 
strength Ogy was set to that value at which the probability of 
failure was 50 percent when a uniaxial load was applied at a 
specified magnitude and direction, as listed in Table II. 
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Ml — Left interface 

• Ml — Top edge 
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- Fiber top edge 
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— Fiber top edge 


Figure 10. — Ml interactive model for the individual sampled points from the unit cell at 
50-percent probability of failure relative to the applied macroscopic composite stresses. 
Also shown are Tsai-Wu and Tsai-Hill models for comparison, (a) Ml without including 
the effect of thermal residual stresses, (b) Ml with the effect of thermal residual stresses 
included. 
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Tsai-Wu 

Tsai-Hill 
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Figure 1 1 . — Ml interactive model for the individual failure modes from the unit cell at 
50-percent probability of failure relative to the applied macroscopic composite stresses. 
Also shown are Tsai-Wu and Tsai-Hill models for comparison, (a) Ml without including 
the effect of thermal residual stresses, (b) Ml with the effect of thermal residual stresses 
included. 
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TABLE VII.— INTERACTIVE MODEL M2 UNIT-SPHERE PARAMETERS 
WITHOUT RESIDUAL THERMAL STRESSES FROM PROCESSING 


Failure mode 

Model and mode 

Weibull 

modulus, 

m v 

Characteristic 

strength, 

aav, 

MPa 

Shetty 

shear- 

sensitivity 

constant, 

C 

Fiber: tensile 3 

L flaw anisotropic 

10.0 

2147 

1.4 

Matrix: transverse tensile 

Isotropic 

10.0 

80.11 

1.4 

Matrix: transverse compressive 

Isotropic 

10.0 

170.6 

NA b 


a Z, longitudinal; one-half angle of anisotropy distribution. A, 1.0°; exponent of sine or cosine 
function, <|>, 0.0. 

'’NA, not applicable. 


TABLE VIII.— INTERACTIVE MODEL M2 UNIT-SPHERE PARAMETERS 
WITH RESIDUAL THERMAL STRESSES FROM PROCESSING 


Failure mode 

Model and mode 

Weibull 

modulus, 

m v 

Characteristic 

strength, 

MPa 

Shetty 

shear-sensitivity 

constant, 

C 

Fiber: tensile 3 

L flaw anisotropic 

10.0 

2128 

1.4 

Matrix: transverse tensile 

Isotropic 

10.0 

63.25 

1.4 

Matrix: transverse compressive 

Isotropic 

10.0 

218.9 

NA b 


a L, longitudinal; one-half angle of anisotropy distribution, A, 1.0°; exponent of sine or cosine 
function, <|>, 0.0. 
b NA, not applicable. 


3 .2.3.1 Model M2 Results 

Figure 12 shows the resultant M2 model failure envelopes 
for the combined failure modes (via Eq. (18)) relative to the 
applied macroscopic composite stresses. As mentioned in 
Section 3.2.3, only three failure modes were used — for fiber 
tensile strength, matrix tensile strength, and matrix compres- 
sive strength — to generate the overall 50-percent probability 
of failure envelopes. Failure envelopes for m v = 10 are shown 
with and without considering the calculated residual stresses 
resulting from processing. In Figure 12 there are several 
conspicuous trends or items. First is that the overall failure 
envelopes show a more rounded appearance (more curvature) 
as opposed to the boxier appearance (lines tending to be more 
linear) of Figure 9. Another is that there was significant 
difference between the envelopes that had and did not have the 
thermal residual stresses included. 

A third important observation is that failure from longitudi- 
nal (axial) tensile loading was predicted to occur at a signifi- 
cantly lower stress than what was measured from the 
experimental rupture data. This was particularly true for the 
failure envelope that included the residual thermal stresses 
from processing. The reason for this was that the matrix was 
predicted to fail (or microcrack) before the fiber bundle failed 
(which was normalized to the longitudinally loaded rupture 
data). If the matrix actually microcracked before fiber failure, 
the composite could be expected to survive intact (with 
reduced stiffness due to the matrix damage) until the load was 


increased sufficiently in order to fail the fibers. Therefore the 
failure envelopes shown in Figure 12 become predominantly 
predictions of matrix cracking and less so predictions of the 
ultimate tensile strength of the composite. Further complica- 
ting the picture in this example was that, depending on the 
multiaxial stress state, the composite may or may not cata- 
strophically fail. For example, under pure uniaxial tensile 
loading transverse to the fibers, the failure should be cata- 
strophic; however, under pure uniaxial tensile loading longitu- 
dinal to the fibers, only matrix failure would occur and not 
necessarily with ultimate catastrophic failure. The difficulty 
then becomes knowing under what combination of transverse 
loading to tensile loading will the composite either catastroph- 
ically fail or remain intact. This indicates the need for progres- 
sive damage modeling, in particular when the composite 
contains laminated plies. 

A fourth important observation in Figure 12 is that failure 
from longitudinal uniaxial compressive loading was predicted 
to occur at much higher levels of compressive stress than the 
experimental results indicated. This was because the fiber 
absorbed much of the loading in the unit cell (the unit-cell 
model did not consider fiber buckling or kinking), and hence 
greater load was required in order for the matrix to reach a 
critical failure stress — a failure stress level that was deter- 
mined (or normalized) from the transverse compressive 
loading on the composite. As discussed in Section 2.3.1, other 
mechanisms such as microbuckling and kinking likely played 
a role in the composite longitudinal compressive-failure 
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Figure 12. — M2 (isotropic matrix strength) interactive model for 
the combined failure modes from the unit cell at 50-percent 
probability of failure relative to the applied macroscopic com- 
posite stresses. Only three failure modes are used — for fiber 
tensile strength, matrix tensile strength, and matrix compres- 
sive strength. The failure envelopes are shown with and 
without considering the predicted residual stresses resulting 
from processing. For the matrix tensile failure mode, Shetty 
shear-sensitivity constant C =1.4. Also shown are adjusted 
stresses whereby multiaxial stresses in quadrants 2 and 3 
are normalized to the longitudinal compressive strength 
experimental data. Also shown are Tsai-Wu and Tsai-Hill 
models for comparison. 

behavior. These mechanisms would cause a shift or redistribu- 
tion of stress such that the matrix would be carrying a higher 
proportion of the compressive load. As the matrix failed 
(crushed) it would result in further fiber kinking and, conse- 
quently, in more stress being transferred to the matrix and so 
on — eventually resulting in composite failure. This implies 
that a multiaxial matrix failure criterion for compression may 
be used to help predict composite failure where a longitudinal 
compressive load is present, provided that fiber kinking (or 
fiber offset) is accounted for in the unit-cell model. As 
previously mentioned, this line of inquiry is beyond the scope 
of this report. However, if the actual magnitude of the matrix 
stresses was higher than what was indicated from the unit-cell 
FE model, a simple adjustment to these stresses such that the 
failure envelope is normalized to the longitudinal failure load 
in compression may produce an acceptable (but a more 
phenomenologically based) failure envelope. This was 
attempted in Figure 12 where the two failure envelopes (with 


and without residual thermal stresses) were adjusted by 
multiplying the longitudinal applied compressive stress for a 
given failure point in the envelope (only in quadrants two and 
three) by the ratio of the experimentally determined strength 
from uniaxial compressive loading in the longitudinal direc- 
tion to that of the predicted value for that same loading. The 
adjusted-stress failure envelopes more closely tracked the 
Tsai-Hill failure criterion than the Tsai-Wu failure criterion — 
except at the highest levels of applied transverse compressive 
stress. If this stress adjustment methodology were actually 
used as a failure criterion for composites, it would have to be 
altered to truncate the failure envelope at the higher levels of 
applied transverse compressive stress to the value of the 
failure stress for pure uniaxial compressive loading in the 
transverse direction (145 MPa in this case), in order to be 
more conservative than the adjusted failure envelope. 

3.2.3.2 For M2: Tracking the Response of Individual 
Sampling Points and Shear Sensitivity 

Finally, it is worth noting that the curvature shown in the 
fourth quadrant better follows the qualitative trend of the 
experimental data than does the prediction in Figure 9 for Ml. 
Tracking of the specific failure modes at the sampled unit-cell 
points is further considered in the next section for model M3 
(and not here for the sake of brevity). However, examining the 
role of shear sensitivity (the C parameter) for the matrix 
tensile-failure mode for M2 (isotropic strength response matrix) 
is worthwhile here (because the further confounding effect of an 
anisotropic strength response matrix is avoided). Referring back 
to the results from Figure 5, the predicted strength response in 
the fourth quadrant was shown to be sensitive to the value of 
C. It was not clear if that same sensitivity in the fourth quad- 
rant also applied in this case. Therefore, Figure 13 was prepared 
showing this effect of parameter C for the 50-percent probabil- 
ity of failure response of the matrix tensile-failure mode for the 
individual sampled points of the matrix from the unit cell. 
Failure models such as Tsai-Wu and Tsai-Hill do not make 
these distinctions regarding location of failure initiation in the 
unit-cell model. Shown are the contour lines for C of 0.82, 1.4, 
and 100 for a highly shear-sensitive, mildly shear-sensitive, and 
shear-insensitive response, respectively. As shown in Figure 5, 
the influence of C on the failure response was smaller for the 
M2 model than for the unit-sphere modeling for graphite. This 
trend was borne out again in the M3 model, which is described 
next. Also in Figure 13, the dominant location (controlling the 
overall failure response in Fig. 12) for the tensile-failure mode 
shifts from the top edge and top interface in the fourth quadrant 
to the left edge and left interface in the first quadrant, and this 
transition seems to occur very near the longitudinal stress axis 
marking the transition between the first and fourth quadrants. 


NASA/TP— 2013-217749 


24 



O Experimental 
— Tsai-Wu 

data 

■ • • Tsai-Hill 

C 

♦ Left edge 

1.4 

X Left interface 

1.4 

• Top edge 

1.4 

Top interface 

1.4 

— • Left edge 

.82 

— • Left interface 

.82 

— Top edge 

.82 

• — Top interface 

.82 

— • Left edge 

100 

■ • • Left interface 

100 

— • Top edge 

100 

■ • • Top interface 

100 


Figure 13. — M2 model of the matrix tensile failure mode versus the applied macroscopic com- 
posite stress showing the effect of the Shetty shear-sensitivity constant, C, on the 50-percent 
probability of failure response of the individual sampled points of the matrix in the unit cell. 
Shown are the contour lines for C of 0.82, 1 .4, and 1 00, for a highly shear-sensitive, a mildly 
shear-sensitive, and shear-insensitive response, respectively. Only the loading case without 
considering residual thermal stresses from processing is shown. Also shown are Tsai-Wu and 
Tsai-Hill models for comparison. 


3.2.4 Model M3: Anisotropic Matrix Material 

In this model variation, designated as M3, the matrix was 
assumed to be mildly anisotropic in strength response. This 
was done specifically so that the matrix did not fail (or 
microcrack) prior to fiber bundle failure under longitudinal 
tensile loading. As with M2, this model only used the unit-cell 
microstresses. Two model variants are described: M3a and 
M3b. M3a used the tensile fiber failure mode described in Ml 
but excluded the compressive fiber failure mode (similar to 
what was done with M2). M3b used the same tensile- and 
compressive-stress fiber failure modes as Ml. In M3a and 
M3b, the matrix tensile-failure mode was assumed to have a 
mildly anisotropic strength response; however, the matrix 
compressive-failure mode used the same isotropic strength 
response model of M2. 

Table IX lists the parameters for the failure modes that do 
not include the effect of the residual thermal stresses from 


processing, and Table X lists the parameters for the failure 
modes that include the effects of residual stresses from 
thermal processing. Only results with the Shetty shear- 
sensitive parameter C= 1.4 for the tensile-failure mode are 
listed. As previously described with Ml, the characteristic 
strength Oq V was set to a value whereby the probability of 
failure was 50 percent when a uniaxial load was applied at the 
magnitude and direction specified in Table II. For the matrix 
tensile-failure mode, the parameters were set for an applied 
tensile load that was parallel to the fiber direction at an 
arbitrarily set value of 1500 MPa (a value not listed in 
Table II), which was sufficiently stronger, for the purpose of 
this example, than the fiber bundle strength for the applied 
macroscopic load of 1280 MPa. The degree of strength 
anisotropy was set so that for an applied uniaxial tensile load 
transverse to the fiber, the strength was 40 MPa (see Table II) 
at a probability of failure of 50 percent. 
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TABLE IX.— INTERACTIVE MODEL M3a AND M3b UNIT-SPHERE PARAMETERS 
WITHOUT RESIDUAL THERMAL STRESSES FROM PROCESSING 


Failure mode 

Model and mode 

Weibull 

modulus, 

m v 

Characteristic 

strength, 

O0F> 

MPa 

Shetty 

shear-sensitivity 

constant, 

C 

Fiber: tensile a-b 

L flaw anisotropic 

10.0 

2147 

1.4 

Fiber: compressive 3 ' 1 ’’ 1 ' 

L flaw anisotropic 

10.0 

1346 

- 

Matrix: transverse tensile d 

T critical strength anisotropic 

10.0 

111.8 

1.4 

Matrix: transverse compressive 3 

Isotropic 

10.0 

170.6 

- 


a The constants y, and r are not used. 

b L, longitudinal; one-half angle of anisotropy distribution. A, 1.0°; exponent of sine or cosine function, (]), 0.0. 
“Failure mode used with model M3b. 

d Set for 1500-MPa longitudinal strength at 50-percent probability of failure; one-half angle of anisotropic 
distribution for critical mode I stress intensity factor, q, 10°; exponent of sine or cosine function, y, 1.0; 
constant, r, 1 .776. 


TABLE X.— INTERACTIVE MODEL M3 A AND M3B UNIT-SPHERE PARAMETERS 
WITH RESIDUAL THERMAL STRESSES FROM PROCESSING 


Failure mode 

Model and mode 

Weibull 

modulus, 

m v 

Characteristic 

strength, 

Oev, 

MPa 

Shetty 

shear-sensitivity 

constant, 

C 

Fiber: tensile 3 ' 15 

L flaw anisotropic 

10.0 

2128 

1.4 

Fiber: compressive 3 ’^ 

L flaw anisotropic 

10.0 

1358 

- 

Matrix: transverse tensile 1 * 

T critical strength anisotropic 

10.0 

127.7 

1.4 

Matrix: transverse compressive 3 

Isotropic 

10.0 

218.9 

- 


a The constants y, and r are not used. 

b L, longitudinal; one-half angle of anisotropy distribution, A, 1.0°; exponent of sine or cosine function, <|>, 0.0. 
“Failure mode used with model M3b. 

d Set for 1500-MPa longitudinal strength at 50-percent probability of failure; one-half angle of anisotropic 
distribution for critical mode I stress intensity factor, 10°; exponent of sine or cosine function, y, 1.0; 
constant, r, 2.959. 


The anisotropic strength response was achieved by using the 
critical- strength and stress-intensity-factor anisotropy model 
described by Equations (16) and (17) for an equatorial-belt 
distribution (see Fig. 3) of angular width i ; T = 0.1745 radians 
(10°) and shape exponent y T = 1.0. The values for parameters 
Z, T and y T were chosen arbitrarily, but they were set so that a 
rather abrupt variation of cti c with orientation occurred. This 
was meant to physically simulate a weaker interface between 
the fiber and the matrix or to be indicative of flaws having a 
lower average strength (or larger average length) normal to (or 
nearly normal to) the fiber axis orientation in comparison to 
other orientations. A simple interval halving search algorithm 
was used to find the values for r T for a user-specified level of 
failure probability and multiaxial stress state. Tables IX and X 
list these values without considering residual thermal stresses 
and considering residual thermal stresses, respectively. 

The amount of strength anisotropy in the matrix is seen by 
comparing the characteristic strength Oer between models M3 
and M2. For example, for the characteristic strength without 


residual stresses, this difference was 111.8 MPa for the 
anisotropic strength matrix versus 80.11 MPa for the isotropic 
matrix (in Tables IX and VII, respectively) for a 1.396 
strength ratio. For the value of Uqv that included the effects of 
residual thermal stresses, this difference was 127.7 MPa for 
the anisotropic strength matrix versus 63.25 MPa for the 
isotropic matrix (in Tables X and VIII, respectively) for a 
higher strength ratio of 2.019. The reason for the increased 
stress ratio can be seen in Figure 12, where a larger amount of 
adjustment is required for the case with the thermal residual 
stresses than for the case without the thermal residual stresses 
in order to normalize the longitudinal tensile strength response 
to the 1500 MPa value. 

3.2.4. 1 Model M3 Results 

Figure 14 shows the resultant failure envelopes for models 
M3a and M3b for the combined failure modes (via Eq. (18)) 
considering and not considering the residual stresses from 
thermal processing. The M3a model is also shown with and 
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Figure 14. — M3a and M3b (anisotropic tensile failure mode for matrix) interactive models for the 
combined failure modes from the unit cell at 50-percent probability of failure relative to the applied 
macroscopic composite stresses. Also shown are adjusted stresses whereby multiaxial stresses in 
quadrants 2 and 3 are normalized to the longitudinal compressive strength experimental data. Model 
M3b is the same as M3a with the addition of the simple longitudinal compressive failure criterion. All 
failure envelopes are shown with and without considering the predicted residual stresses that resulted 
from processing. For the matrix tensile failure mode, Shetty shear-sensitivity constant C = 1 .4 was 
assumed. Also shown are Tsai-Wu and Tsai-Hill models for comparison. 


without the adjusted stresses (as per M2). From the figure it is 
observed that allowing the strength response of the matrix in 
tension to be anisotropic brought the failure envelopes in 
better alignment with the experimental data in the fourth 
quadrant (both with and without including the effects of 
residual thermal stresses) and closer to (or between) the results 
of the Tsai-Wu and Tsai-Hill failure criteria. The M3 and M2 
models illustrate the caution that must be exercised when 
applying Equation (18): it is important to know what the 
predominant failure modes are and whether or not they imply 
ultimate composite failure. 

In quadrant one of Figure 14, the results from the M3 model 
that included and did not include residual stresses were in 
closer agreement than were the results from the Ml model in 
Figure 9, although the M3 model results were also slightly 
more conservative in comparison to the experimental data 
shown in that quadrant. The predicted envelopes fell between 
the Tsai-Wu and Tsai-Hill failure envelopes in that quadrant. 

There are also some noteworthy features in quadrants two 
and three. With the M3a model (with and without the effect of 
residual stresses), the failure envelopes for the unadjusted 
stresses do show some difference in shape relative to the 
failure envelope of the M2 model in Figure 12. This was 
particularly surprising in quadrant three because both the M2 
and M3a models used the same isotropic strength response 
model for the matrix in compression. This discrepancy is 
explored further when the individual failure modes are 
examined. The difference in shape was also transferred to the 


adjusted stresses, where the predicted results were somewhat 
more nonconservative in comparison to the Tsai-Hill criterion 
than were the results in Figure 12. Nevertheless, the trend 
remained that the adjusted stresses tracked much closer to the 
Tsai-Hill criterion than to the Tsai-Wu criterion. Quadrants 
two and three of Figure 14 also show the effect of the fiber 
failure mode in compression for model M3b. Although 
employing this failure mode was a somewhat trivial exercise, 
it was useful to compare the predictions of this failure mode to 
the predictions with the adjusted stresses. Overall, Figure 14 
showed that the M3a model with the stresses adjusted in 
quadrants two and three followed the experimental data 
reasonably well and also that the failure envelopes (with and 
without considering residual thermal stresses) tracked rela- 
tively close to the Tsai-Hill failure criterion. 

3.2.4.2 For M3: Tracking the Response of Individual 
Sampling Points and Failure Modes 

Figure 15 shows the matrix tensile and compressive-failure 
modes for the individual sampled points in the unit cell for the 
M3a model without consideration of the effect of residual 
thermal stresses from processing. This figure reveals that the 
dominant failure modes were relative to the stress quadrant as 
well as the sampling points that were the source for this 
behavior. To see the contribution of the fiber failure modes, 
refer back to Figures 10 and 11. Some interesting trends are 
revealed in Figure 15. For example, the top edge sampling- 
point is shown to contribute significantly to matrix failure in 
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Figure 1 5. — M3a model of individual points for 50-percent probability of failure for the matrix 
compressive (C) and matrix tensile (T) failure modes without consideration of the effect of 
the residual stresses from thermal processing. For the matrix tensile failure mode, Shetty 
shear-sensitivity constant C = 1.4 was assumed. Also shown are Tsai-Wu and Tsai-Hill 
models for comparison. 


the third quadrant (also compare this with Fig. 10) in both the 
tensile and compressive-failure modes. Also, in the third 
quadrant for the most extreme compressive loading case (the 
combination of the highest longitudinal and transverse 
compressive stresses acting simultaneously), a tensile-failure 
mode is indicated. 

In the fourth quadrant of Figure 14, the failure envelope 
showed a curved path (or rounded corner). Figure 15 indicates 
that a major factor in this behavior was the compressive- 
failure mode at the left-edge sampling point. This contrasts 
with Figure 5 for graphite where the tensile-failure mode 
played a more prominent role in this (fourth) quadrant as 
evidenced by the role that the shear-sensitivity parameter C 
played in the resultant shape of the failure envelope. Another 
interesting observation in the fourth quadrant was that the 
tensile-failure mode at the top edge and top interface points 
played a significant role in matrix failure at the higher levels 
of applied longitudinal tensile stress. In the M2 model shown 
in Figure 12, these failure contour lines shifted to the left 
(became weaker) and resulted in the matrix being predicted to 
microcrack before fiber bundle failure would likely occur. 
Also of note in Figure 15 was that the location for the tensile- 
failure mode shifts from the top edge and top interface in the 
fourth quadrant to the left edge and left interface in the first 
quadrant and that this transition seems to occur very near the 
portion of the longitudinal stress axis that marks the transition 
between the first and fourth quadrants. This was also indicated 
to occur in the M2 model in Figure 13. Also, in the first 
quadrant of Figure 15 it can be seen that the prediction of 


matrix failure correlated well with the experimental data. The 
greater conservatism seen in the failure envelope in Figure 14 
was, therefore, due to the added effect of the fiber tensile- 
failure mode on the overall failure probability via 
Equation (18). The experimental data in the first quadrant 
were actually an adjustment made by Al-Khalil et al. (1996) 
that added a transverse tensile-stress component to what was 
originally intended as pure longitudinal tensile loading. A 
fiber failure mode was indicated for these data. However, it is 
also interesting to note that Al-Khalil et al. (1996) initially had 
difficulty suppressing catastrophic matrix failure transverse to 
the fiber under longitudinal tensile loading — perhaps illustrat- 
ing some sensitivity to the transverse loading component. 

Figure 16 shows the effect of parameter C for the 
50-percent probability-of-failure response of the matrix 
tensile-failure mode for the individual sampled points of the 
matrix from the unit cell. The individual points are shown 
(as opposed to the overall failure probability) so that the rela- 
tive contribution and sensitivity to C of each point location 
on the overall failure probability can also be seen. Shown are 
the contour lines for C of 0.82, 1.4, and 100 for a highly 
shear-sensitive, a mildly shear-sensitive, and a shear- 
insensitive response, respectively. Figure 16 shows that the 
influence of C on the failure response was small for the M3 
models in comparison to the unit-sphere modeling for graphite 
shown in Figure 5. This confirmed that adding strength 
anisotropy to the M3 model did not appreciably affect the 
behavior of C in comparison to the results from M2 shown in 
Figure 13 for the isotropic strength response matrix. 
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Figure 16. — M3a model of matrix tensile failure mode showing the effect of the shear sensitivity 
parameter C (without thermal residual stresses). Also shown are Tsai-Wu and Tsai-Hill 
models for comparison. 


Overall, the M3 models provided better correlation to the 
experimental data than did the Ml and M2 unit-sphere models 
(depending on when matrix microcracking initiated). Better 
correlation to the experimental data in the fourth quadrant 
could have been achieved if the strength anisotropy in the 
matrix had been reduced by making it weaker in the fiber 
direction. However, doing so would have implied that matrix 
failure would become more predominant over the fiber failure 
mode at the highest levels of applied longitudinal stress, and 
this was not pursued. 

4.0 Summary and Conclusions 

This report introduced the unit-sphere methodology to the 
problem of multiaxial strength response for brittle or quasi- 
brittle anisotropic and composite materials. This required the 
development of modeling extensions that allowed for (1) flaw- 
orientation anisotropy, whereby a preexisting microcrack has a 
higher likelihood of being oriented in one direction over 
another direction; and (2) critical strength or fracture tough- 
ness anisotropy, whereby the level of critical strength <3\ c or 
fracture toughness K lc for mode I crack propagation changes 
with regard to the orientation of the microstructure. This was 
done for transversely isotropic strength response materials for 
the tensile mode of failure. In addition, a simple shear-stress- 
based unit-sphere failure criterion was introduced to account 
for the compressive mode of failure for isotropic-strength- 
response materials. A central principle to this model was that 
various failure modes controlled the material strength response 


and that these modes are stochastic and described by a 
Weibull distribution. The probability of material survival 
(the inverse of probability of failure) was the product of the 
probabilities of survival for the various failure modes. This 
assumed that weakest-link behavior was operative with no 
interactivity between failure modes. 

The unit-sphere model for anisotropic strength response was 
applied to nuclear-grade IG-110 and H-451 graphite and to a 
unidirectional-fiber-reinforced polymer-matrix composite 
under biaxial loading situations (uniaxial loads applied 
orthogonal to one another). First, this was done to demonstrate 
the methodology for a macroscopically homogeneous material 
and, subsequently, for a heterogeneous material system. 

For the nuclear-grade graphite, failure envelopes for con- 
stant levels of probability of failure for biaxial loading were 
generated for an isotropic strength response material (the 
IG-110 grade) and a mildly anisotropic strength response 
material (the H-451 grade). The experimental data for pure 
uniaxial loading were used to normalize or calibrate the 
models. The Shetty shear-sensitive mixed-mode fracture 
criterion with an assumed penny-shaped crack modeled the 
flaws in the materials for the tensile-failure mode, and the sim- 
ple shear-based failure criterion was used to predict the 
response for the compressive-failure mode. Good correlation 
was achieved for the H-45 1 rupture data, which was available 
only for the first and fourth biaxial-loading quadrants (the 
isotropic strength response for the near-isotropic IG-110 grade 
material was shown mainly for comparison purposes). Failure 
envelopes that assumed highly shear-sensitive (to mode II) 
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and shear-insensitive flaws were also shown. The model 
indicated that the level of shear sensitivity significantly 
affected the predicted strength response in the fourth loading 
quadrant (for applied tensile and compressive stresses) and 
that a mildly shear-sensitive flaw response correlated best to 
the experimental data. 

For the unidirectional composite problem, the test case from 
the WWFE for the unidirectional lamina under combined 
longitudinal and transverse loading was chosen. An FE model 
of a continuous-fiber-reinforced composite unit cell with 
square packing was examined. The thermal residual stresses 
that develop from processing (during cooldown) were also 
accounted for in the linear-elastic stress analysis. The stress 
analysis revealed that four locations (two of the locations were 
coincident with the fiber matrix interface) contained the 
highest magnitude stresses. These four locations were subse- 
quently used in the context of a maximum-stressed-point 
failure criterion to predict the probability of failure of the 
composite. Failure envelope predictions (at 50-percent 
probability of failure) from the unit-sphere model were 
compared with the experimental results of the WWFE test 
case as well as to the well-known Tsai-Wu and Tsai-Hill 
failure criteria. 

Four failure modes were assumed to be active for the 
composite material: (1) fiber failure in tension, (2) fiber failure 
in compression, (3) matrix failure in tension, and (4) matrix 
failure in compression. Three variations of the anisotropic unit- 
sphere model were examined — designated as Ml (flaw and 
failure plane anisotropy), M2 (isotropic matrix strength), and 
M3 (anisotropic matrix strength). The major assumption behind 
the Ml model was that the flaw planes in the matrix were either 
oriented parallel to the fiber axis or that strength was controlled 
by global failure (action) planes that ran parallel to the fiber 
axis. A noninteractive version of Ml that used only applied 
global macroscopic stresses also was investigated. The resulting 
failure envelopes (with thermal residual stresses and without 
thermal residual stresses) did not agree with the experimental 
data well in the fourth stress quadrant. It was shown that a 
transverse stress component in the matrix that became more 
compressive as longitudinal load increased affected the 
compressive-failure mode and was responsible for the shape of 
the failure envelope in the fourth quadrant. The noninteractive 
Ml version corresponded closely to a stochastic version of the 
maximum principal-stress failure criterion. 

Conversely, the M2 model assumed an isotropic strength 
matrix, and the compressive fiber failure mode was not 
considered in the overall failure probability calculation. The 
resultant failure envelopes took a distinctly more rounded 


shape than the results from Ml because of the added contribu- 
tion of the longitudinal stress component in the matrix. Major 
findings indicated that (1) the matrix was predicted to micro- 
crack before the fiber bundle failed, (2) matrix compressive 
strength was much greater than the composite longitudinal 
compressive strength, and (3) there were significant differ- 
ences in the predicted strength response when residual thermal 
stresses were included as opposed to when they were not. The 
predicted failure envelopes did not correlate well with either 
the experimental data or the Tsai-Wu or Tsai-Hill failure 
criteria. Normalizing (pinning) the longitudinal uniaxial 
compressive response to the experimental data enabled failure 
envelopes in biaxial stress quadrants two and three to compare 
more favorably with the Tsai-Hill failure criterion results in 
these two quadrants. 

The M3 model assumed that the matrix had a mildly aniso- 
tropic strength response for the tensile-failure mode so that 
matrix microcracking would be predicted to occur at an 
applied longitudinal tensile stress that was higher than the 
experimental strength of the fiber bundle. For M3 the overall 
failure envelopes tracked quite well to the limited amount of 
experimental data that were available and clearly provided 
better correlation than the Ml and M2 models. 

The development of the anisotropic unit-sphere methodol- 
ogy was an attempt to provide an improved mechanistic basis 
to the problem of predicting composite strength under multi- 
axial loading compared with phenomenological-based poly- 
nomial formulations such as Tsai-Wu, Tsai-Hill, and Flashin, 
among others. This model is not as simple and straightforward 
to use as these other failure criteria; however, the value of this 
work is in the insights that are provided about stress-state 
interactions, possible failure modes, and the physical assump- 
tions behind predicted behaviors. 

Planned future work will incorporate unit-sphere methodol- 
ogy into the NASA Glenn Research Center micromechanics 
analysis code/generalized method of cells (MAC/GMC) 
analysis tool. This will involve the full-field stress solution of 
the unit cell, including stressed-volume size effects. The ever- 
expanding capabilities of MAC/GMC combined with the time- 
stepping (and time- and cycle-dependent) transient reliability 
analysis capabilities of the unit-sphere algorithm conceivably 
will allow for stochastic-based progressive damage modeling 
capability that takes into account time- and cycle-dependent 
material degradation in a computationally efficient manner. 

Glenn Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, August 1, 2013 
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Appendix — Symbols 


area 

critical crack length or crack radius 
Shetty shear-sensitivity coefficient 
constant 

Young’s modulus of elasticity 
normalized anisotropy function of K k: or 
ct Ic as a function of angles a and P 
Heaviside step function 
mode I stress-intensity factor 
critical mode I stress-intensity factor 

Batdorf uniaxial stress-state normality 
constant for the volume-flaw failure 
mode (also known as the Batdorf nor- 
malized crack-density coefficient) 

Weibull modulus 

Weibull modulus for the volume-flaw 
failure mode for compressive stress 
states 

Weibull modulus for the volume-flaw 
failure mode 

number of failure modes 
probability of failure (P/= 1 —P s ) 
probability of failure of material volume 
probability of failure of a crack with a 
strength between Gi c and (gi c + Acti c ) 
in AF 

reliability or probability of survival 
(P s =l-P f ) 

probability of existence of a crack with 
strength between CTi c and (gi c + Acti c ) in 
an incremental volume 

probability that a crack of critical 
strength will be oriented in a particular 
direction such that it will grow and 
cause failure 

probability that a crack of critical 
strength 0 \ eqc is oriented in the range 
between a and (a + da) and between p 
and (P + dP) 

constant (ratio) or parameter in K\ c 
anisotropy function for polar-cap 
distribution 

constant (ratio) or parameter in K lc 
anisotropy function for equatorial-belt 
distribution 


V 
V e 

x, y, z 

Y 

a, P 


a, 

Y 


lL 


Y T 


m 

C(a, P) 


n [( g i(Y/c) 


0 

A 


A L 


V 

% 


volume 

effective volume 

location in the body of the structure; 
Cartesian coordinates 

crack-shape geometry factor 

orientation angles or angular coordinates, 
where the direction normal to the plane of 
the microcrack is specified by the radial 
line defined by a and P in stress space 
thermal expansion coefficient 

constant or parameter in K lr anisotropy 
function representing the exponent of 
the sine or cosine function 

constant or parameter in K lc anisotropy 
function representing the longitudinal 
(polar-cap) distribution 

constant or parameter in K lc anisotropy 
function representing the transverse 
(equatorial-belt) distribution 

C, as a function of angle a describing the 
anisotropy of flaw orientation 

function describing the anisotropy of 
the flaw orientation where the normal 
direction to the flaw plane is given by 
angles a and P 

crack-density function for equivalent 
mode I strength, < 5 \ eqc , of a flaw: number 
of flaws per unit volume with strength 
equal to or less than Gi e?c 

angle measured from the longitudinal 
loading axis 

constant or parameter in flaw- 
orientation anisotropy function repre- 
senting one-half of the total angular 
extent of the anisotropy distribution 

constant or parameter in flaw-orientation 
anisotropy function representing the 
longitudinal (polar-cap) distribution 

constant or parameter in flaw-orientation 
anisotropy function or in K lc anisotropy 
function representing the transverse 
(equatorial-belt) distribution 

Poisson’s ratio 

constant or parameter in K lc anisotropy 
function representing one-half of the 
total angular extent of the anisotropy 
distribution 
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constant or parameter in AT| ( anisotropy 
function representing the longitudinal 
(polar-cap) distribution 

4*7- 

constant or parameter in flaw-orientation 
anisotropy function representing the 
transverse (equatorial-belt) distribution 

\ T 

constant or parameter in K u anisotropy 
function representing the transverse 
(equatorial-belt) distribution 

£2(X, G J eqc ) 

area of a solid angle projected onto a 
unit radius sphere in three-dimensional 
stress space for which c Ie? > <3\ eqc from 

x 

summation function; applied far-field 
multiaxial stress state 


an applied multiaxial stress state X 

CT 

applied uniaxial stress 

Superscript 


Gie 

®Ic,x» ^Ic,y’ °Ic,z 

critical mode I strength 

orthogonal critical strength components 

normalized by G Ic>max 

Subscripts 

normalized 

^ 1 (Y/ 

equivalent, or effective, stress 

e 

compression 

^Ieg.max 

maximum value of a\ eq over the unit 

c 

critical 

sphere from the applied multiaxial 

eq 

equivalent, or effective 


stress X 

1 

mode I 

G | ( Y/C 

critical equivalent, or effective, strength 

i 

/ th value or I th term 

ct^x, y, z, a, P) 

G n 

critical equivalent, or effective, strength 
of a flaw located at coordinates x, y, 
and z and oriented at angles a and (1 
applied far-field stress component 

L 

max 

T 

longitudinal 

maximum 

transverse 

GoV 

normal to a crack face 

Weibull scale parameter for the volume- 
flaw failure mode normalized to unit 
volume 

V 

Definitions 

volume or a volume -based property 
(e.g., indicates volume-flaw analysis) 

°x, t7y, °z 

orthogonal stress components expressed 
relative to a global coordinate system 

fast fracture 

component rupture in the absence 
of slow crack growth where 


Weibull characteristic strength (volume- 
flaw failure mode) 


strength is strictly controlled by 
the fracture toughness and the 
size, distribution, and orientation 


Weibull characteristic compressive 
strength (volume-flaw failure mode) 

j 

of inherent flaws 
longitudinal 

T 

shear stress acting on the oblique plane 

L, 


whose normal is determined by angles 

mode 1 

crack-opening mode 


a and (3, which represents the applied 
far-field shear stress on a crack face 

mode 11 

crack-sliding mode (in-plane 
shear) 

"txy, "tyz> Tzx 

shear stress components expressed 
relative to a global coordinate system 

mode 111 

crack-tearing mode (out-of-plane 
shear) 

<l> 

constant or parameter in flaw- 

T 

transverse 


orientation anisotropy function repre- 
senting the exponent of the sine or 
cosine function 

transient reliability analysis predicting the probability of 

survival of a component while 
accounting for loads and temper- 

<l>i 

constant or parameter in flaw-orientation 
anisotropy function representing the 
longitudinal (polar-cap) distribution 


atures that can vary over time 
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